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Space-time double correlations and spectra in a 
turbulent boundary laver 


By A. J. FAVRE, J. J. GAVIGLIO, and R. DUMAS 
Laboratoire de Mécanique de l’ Atmosphére, Université de Marseille 


(Received 3 September 1956 and in revised form 30 January 1957) 


SUMMARY 

This paper describes the results of an experimental investigation 
of the turbulent boundary layer on a flat plate with zero pressure 
gradient, carried out at the Laboratoire de Mécanique de 
l’ Atmosphere (1.M.F. Marseille). ‘Transition to turbulent flow 
was obtained either by increasing the preturbulence upstream of 
the plate by means of a grid, or by a series of emery paper 
roughness elements beginning at the leading edge. ‘The measure- 
ments were of the space-time double correlations, i.e. double 
velocity correlations with both spatial separation and time delay, 
from which isocorrelation lines for optimum delay can be drawn, 
and of the spectrum function, all for the longitudinal component 
of the velocity fluctuation. 


1. EXPERIMENTAL APPLIANCES 

Most of these have been described in previous papers (Favre, Gaviglio 
& Dumas, 1952, 1953 b, 1954a, 1954b). 

The wind-tunnel has a complete return circuit, and is operated by two 
contra-rotating fans whose pitch is adjustable during operation. ‘The fans 
are driven by two direct-current motors, fed by a 60 H.P. regulated ignitron 
rectifier, making it possible to change the rotation speed. ‘The test section 
is 80cm x 80cmx270cm. ‘There is a smooth contraction with an area 
ratio of 15; it is preceded by a settling chamber, in which five damping 
screens are set, and by a rapid expansion diffuser, with a screen analogous 
to that of the low-turbulence wind-tunnel of the R.A.E. (Squire, Winter & 
Barnes 1947). ‘The natural intensity ratio of the turbulence in the test 
section is 0-00038 for a speed of 12-20 m/sec. 

A flat plate made of glass (Favre, Gaviglio & Dumas 1954a), 210 cm 
in length and 80 cm in width, and with a profiled leading edge, was hung 
in the test section. ‘The measurements have been made, as previously, 
under the lower side of the plate at distances from the leading edge of 
about 80 cm and 194 cm. 

y’ and y are distances from the walls of the wires A (fixed) and B 
(adjustable); Y, = y’—y (figure 1); CB =X, is their separation in the 
direction of the general stream velocity; their separation Y, parallel to the 
plate and normal to the stream is zero. 7 is the interval of time added to 
the instant of the signal corresponding to the fixed wire A. 
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‘The hot-wire probes used were either of an old pattern (diameter ot 
the body 10 mm—-see figure 1) or of a new type (diameter of the body 8 mm, 
with holders thinned and lengthened); the holders of the upstream probe 
were made lozenge-shaped in order to lessen the effects of the wake. 

The platinum hot wire was made by the wollaston process with an 
etched length 7, in diameter and about 1 mm in length. In the position 
z = 194 cm, the new type of probe holder was used. ‘This holder, fixed to 
the upper side of the plate, passed round the trailing edge and was supported 
on the lower side. It allowed variation of X,, Xo, X3. 


jemery cloth n°00 
lemery {cloth n°2 
cloth n°5 
Figure 1. Experimental arrangement: transition by roughness on leading edge, 
or by preturbulence downstream of a grid. 


‘Transition to turbulent flow in the boundary layer was obtained either 
by increasing the preturbulence of the stream by means of a biplane grid 
of mesh M/ = 1 in. with rods of circular section (d = 0-5 cm) (Favre, Gaviglio 
& Dumas 1952) set at 20cm upstream of the leading edge of the plate 
(figure 1), or by inserting, starting at this leading edge (figure 1), roughness 
in the form of emery cloth (no. 00 up to s = 20cm, then no, 2 to zs =30em, 
and then no. 5 to s = 40 cm). 

The principal measuring appliances were the hot-wire anemometers 
reproducing the oscillations, within +4°, in energy from 1 to 3500 hertz, 
the spectral analyser with a selectivity of 45°, of the median frequency 
and a measuring band extending from 1 to 5000 hertz, and the apparatus 
for measurement of time and space correlations (Favre 1946, 1948; Favre, 
Gaviglio & Dumas 1953 a) with a band-pass extending from 1 to 2500 hertz. 
‘This latter apparatus records simultaneously, and reproduces with a relative 
delay, the two voltages coming from the two anemometers, and measures 
the time-correlation coefficient. 


2. CHECKS ON THE MEASUREMENTS 


We tried to check the influence of the disturbances due to the upstream 
wire (platinum wire, silver coating, needles, body of the probe) on the 
precision of measurements of the intensities of turbulence, the local mean 
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velocity and the space-time correlation (Favre, Gaviglio & Dumas 1955). 
With the old type of probes the disturbances in unfavourable conditions 
(upstream wire near the wall) may have reached 10°, of the turbulence 
intensities, 5°,, of the mean velocity, and 0-05° , of the space-time correlation. 
With the new type probes these disturbances seemed to be reduced to 4°, 
for intensities and to 1°, for the mean velocity. 

We have verified, thanks to a method proposed by H. W. Liepmann 
(1952), that the length of the hot wires had a negligible influence on the 
measurements of space-time correlation (for these experimental investiga- 
tions, the minimum value of the dissipation length A was 2:5 mm). 


3. MEAN VELOCITIES 
‘The velocity Vax outside the boundary layer was 12-00 m/sec. In the 
case of the measurements made at s = 79 cm, with hot-wire holders, the 
mean gradient of static pressure is negative and less than 4°% per metre. 
For the measurements made at x = 194cm, with the new holders, the 
mean gradient is negative and less than 1°, per metre downstream of the 
roughness. 
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Figure 2. Mean velocity in the boundary layer on a flat plate downstream of a grid. 


In the case of the boundary layer on a flat plate downstream of a 
turbulence grid, the relative velocities V,/Vmax were determined as a 
function of the distance y from the wall, for x = 12:5, 23-7, 35-5, 48-0, 
66-0, 74-0, 80-3 cm, and the results are shown in figure 2. ‘The thickness 
of the boundary layer at 79 cm from the leading edge was 16°38 mm. 

In the case of the boundary layer on a flat plate with a rough leading 
edge, the velocity profiles were measured, in the absence of the hot-wire 
holders, at zs = 58-5, 79-0, 167-5 and 194-0 cm, and are given in figure 3. 
A supplementary curve has been drawn to show the effect of the presence 
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of the holders. ‘Vhe velocity V’, in the vicinity of the wall was determined 
at s = 194m by means of a hot wire. ‘The thickness of boundary layer 
in this case was 17-5 mm at s = 79 cm, and 34 mm at z = 194 cm. 

In figure + we have represented the displacement thickness 5* and the 
momentum thickness 4 as functions of the distance = from the leading edge. 
‘The roughness set on the leading edge brings about a marked increase in 
the thickness of the boundary layer. 


“4 4 y mm 


Figure 3. Mean velocity in the boundary layer on a flat plate with a rough leading 


edge. 
8 4S%mm 
6 
rough leading edge | 
4 
| 
2 — 
0 40 80 


Figure +. Displacement thickness 6* and momentum thickness 6. 


Figure 5 shows Clauser’s universal velocity profile (Clauser 1954). 
In this figure, Vy = 4 (7) p) with 7, = pI=,, dé ds, and the value of 5 
corresponds to TV, Vmax = 1. In the case of transition obtained by 
preturbulence (in the presence of the old type of wire holders, the mean 
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pressure gradient being 4°, per metre) the data are not wholly consistent 
with the curve, but in the case of transition obtained by roughness (with 
the new type of holders, the mean pressure gradient being 1% per metre) 
the agreement between the data and the curve is fairly good. 
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Figure 5. Comparison of measured velocities with Clauser’s universal profile. 


4. INTENSITIES AND SPECTRA 

In the case of transition by preturbulence, the intensities as functions 
of the distance v6 have been given previously (Favre, Gaviglio & Dumas 
1954a). ‘The intensity (77), reaches 15°, for = 0-06, which may 
be compared with the value 1-8°,, for the intensity of free turbulence at 
a great distance from the wall. In the case of the transition by roughness, 
two series of measurements have been made, at z = 79 cm and = = 194 cm. 
The intensity, as a fraction of the local mean velocity, seems to decrease 
more rapidly at s = 194cm than at s = 79cm when y/6 increases. At 
> = 79 em, it reaches 18", for y 6 = 0-06, that of the free turbulence being 
of the order of 0-045°,.. 

The spectral functions of turbulence energy F(n) (Taylor 1938) are 
obtained as functions of the frequency ”; directly by selective amplification, 
and indirectly by Fourier transformation of the auto correlation curves. 

‘The areas of the curves F() have been chosen so as to be equal for the 
common intervals of frequencies. 

In the case of transition by preturbulence, the spectra have been 
published previously (Favre, Gaviglio & Dumas 1954a). ‘They correspond 
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tos =79cm. ‘The difference between these spectra inside and outside the 
boundary layer is small, except for low frequencies (less than 40 hertz) 
where they change markedly as functions of the distance from the wall, 
thus disclosing an increase of energy in the boundary layer. 

In the case of transition by roughness, unlike the previous case, the 
intensity of turbulence is very low outside the boundary layer (0-045",,). 
Figure 6 shows the spectra measured at 79 cm from the leading edge for 
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Figure 6. Turbulence spectra in the boundary layer on a flat plate with rough leading 
edge. (Rs 14 500, 79 cm). 


different distances y from the plate. Inside the boundary layer the 
evolution at frequencies greater than 80 hertz is slight; for the lower 
frequencies the curves show an increase of energy as v 6 decreases. Outside 
the conventional boundary layer, at v6 = 1-31 and v6 = 1-71 (the factor 
of intermittance there being markedly smaller than 0-04 according to 
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Corrsin & Kistler (1954)), the spectra become different; each of them 
presents a marked maximum. At the distance s = 194 cm from the leading 
edge (figure 7) the results are qualitatively the same, but with lower 
frequencies, which corresponds to the increase of the scale of turbulence 
with the thickness of the boundary layer. Outside the boundary layer, 
for instance, at y/S = 1-35 and y/6 = 1-77, a marked maximum was found 
as in the previous case. 
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with rough leading 


edge. (R; = 28 000, z = 194 cm). 

The intensities of these fluctuations outside the boundary layer are 
0-47°% and 0-17% respectively for y/8 = 1:31 and 1-71 (z = 79 cm), and 
0:38% and 0-13%, for y/5 = 1:35 and 1-77 (x = 194cm). Thus, they do 
not seem to be due to the preturbulence, the intensity of which is of 0-045°%. 
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When we inserted a small thin plate (4-2 cm = 12-5 cm) parallel to the 
wall, at y 6 = 1-25, between the hot wire placed at v/6 = 1-8 (z = 79 cm) 
and this wall, the intensity was reduced by more than half and this 
phenomenon was no longer detected. 

It may be that these fluctuations are induced by random motion at the 
edge of the boundary layer (Phillips 1954). 
ny/Fee ' 


Figure 8. Comparison of turbulence spectra in several boundary layers on flat plates. 


We have compared (figure 8) the spectra obtained for s = 79 cm and 
z= 194 cm 6 having respectively the neighbouring values 0-69 and 0-71), 
with Klebanotf & Diehl’s (1952) spectra. The ordinates and abscissae in 
figure S are F(n) F(e) and nF(e), F(e) being the average value corresponding 
to the lowest frequencies. ‘The results are fairly consistent. 


5. MEASUREMENTS OF SPACE-TIME CORRELATIONS 
First consider the measurements of transverse space-time correlations 
(7, X, = 0, = 0,.X5). 
In the case of transition by preturbulence, the transverse space-time 
correlation R,,(7, X, = 0, X, = 6, X3) is represented in figures 9 to 12. 
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figure 9, y'/d = 0-06, y’ = 1mm, s = 79cm; 
figure 10, vy’ 8 = 0-06, y = Imm, s= 74cm; 
figure 11, y’/5 = 0-24, y = 4mn, == 79 cm; 
figure 12, y’ /d = 0-71, y’ = 12 mn, z= 79 cm. 
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Figure 9. ‘Transverse space-time correlations in the boundary 
downstream of a grid. 
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Figure 10. Transverse space-time correlations in the boundary layer on a flat plate 
16-0 mm). 


downstream of a grid. 
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74 cm, 0 
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One finds that the correlation reaches a maximum for a value 7; of the delay 
which is a function of Y,. Figure 13 summarizes the values of 7;. Note 
that for the distances from the leading edge s =74 cm and z= 79 cm, 
which differ only by 36, the values of 7; are very little different. Note also 
that 7; changes sign with XY, = y’—y. 
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Figure 11. Transverse space-time correlations in the boundary layer on a flat plate 
downstream of a grid. (v’ = 4mm, z = 79 cm, 6 = 16:8 mm). 


In the case of the transition by roughness, the measurements have been 
made at two different stations, at distances z = 79 cm and z = 194 cm from 
the leading edge. ‘The measurements of R,,(7, X, = 0, X, =0, X53) at 
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the distance s = 79 cm are represented in figures 14 to 16, for which the 


data are: 


figure 14, y'/d = 0-06, y = I] mm; 
figure 15, y'/6 = 0-23, y +mm; 


ll 


tigure 16, v6 = 0-69, y’ = 12 mm. 
R 
1X3 =25 4.0.24 
X3 = 50 
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Figure 12. Transverse space-time correlations in the boundary layer on a flat plate, 
downstream of a grid. (\” 12 mm, z 79 cm, 0 16-8 mm). 


The measurements at the distance s = 194 cm are represented in figures 17 


to 19, for which the data are: : 
tigure 17, vy’ 6 = 0-06, y= Ze; 
tigure 1S, y’/d = 0-23, y = Smm; 
figure 19, v6 = 0-69, = 24+ mm. 
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Figure 13. Optimum delays 7 7; for transverse space-time correlations in the 


boundary layer on a flat plate downstream of a grid. 
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Figure 14. ‘lransverse space-time correlations in the boundary layer on a flat’ plate 
with rough leading edge (\’ 


imm 3: 79 cm, 6 17-5 mm). 
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As in the case of transition by preturbulence, one finds that the correlation 
reaches a maximum for a value 7, of the delay which is a function of 3. 
Figure 20 summarizes the values of 7; which change sign with X3. 

The values of 7; obtained at the distances : = 74cm and s = 79cm 
from the leading edge and for a thickness of the boundary layer of the same 
order are only slightly different, whether the transition is due to pre- 
turbulence (figure 13) or roughness (figure 20). Also, the curves for the 
distribution of the corresponding mean velocities are also only slightly 
different. 
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Figure 15. Transverse space-time correlations in the boundary layer on a flat plate 
with rough leading edge (v’ — 4mm, z = 79 cm, 6 = 17:5 mm). 


The measurements of 7; made at the two distances s = 79cm and 
* = 194 cm (distances which differ by about 45 times the average thickness 
of the boundary layer) show that, Y, being fixed, the delay 7; generally 
increases with the distance s from the leading edge, in the direction of 
increasing thickness of the boundary layer. Moreover, this increase of 7; 
is much smaller than the differences of the delays 7, corresponding to the 
transport due to the mean motion, if one considers this transport over the 
distance separating the position s = 79 cm from the position == 194 cm 


(this delay, along the isochronous lines is 7,-—7, = (3—3)(1/V’—1/V)). 
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Figure 16. Transverse space-time correlations in the boundary layer on a flat plate 
- 17-5 mm). 
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‘Table 1 gives this comparison, in the case of the distances z mentioned 
above. 

Measurements of R),(7, X, = 0, X, = 0, X3) have also been made by 
means of a double analyser which admitted only the waves of frequencies 
of about 90 hertz. Whenz =79 cm, v’/5 = 0-69, X4/5 = 0-40 (nF(e) = 0-32), 
we see from figure 16 that the maximum correlation 0-33 is reached for a 
delay 7; of about 0-95 ms. When zs = 194cm, v’/5 = 0-71, X,/5 = 0-38 
(nF(e) = 0-60), figure 19 shows that the maximum correlation 0-22 is 
reached for a delay 7; of about 2:2 ms. So, in both cases, the value of 7; 
is not practically modified, and that of the correlation is increased by 
filtering, but this is not very significant because these frequencies are in 
the middle of the spectrum. 
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Figure 17. Transverse space-time correlations in the boundary layer on a flat plate 
with rough leading edge (v’ = 2mm, z = 194cm, 6 = 34mm). 


We turn now the longitudinal space-time correlations 
= 0, X3 = 0). Let us recall that we tried (Favre, Gaviglio 
& Dumas 1953b) to obtain a verification of Taylor’s hypothesis for the 
turbulence downstream of a grid, and found a satisfactory concordance 
of the autocorrelation curves with the curves of longitudinal space 
correlation, established up to a distance X,/M of about 8, with the 
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correspondence 7’ = Y,/V’; but for distances Y,/.W > 1-5, the correlations 
being smaller than 0-1, the comparison is no longer conclusive. 

Further verification, including measurements of the longitudinal 
space-time correlations (up to.Y,/M = 8) with optimum delay to compensate 
for the transport due to the mean flow, has shown that these correlations 
retain high values, of the order of 0-5, up to great distances Y, W = 5 
to 8, these values not being, however, equal to unity. 
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Figure 18. ‘Transverse space-time correlations in the boundary layer on a flat plate 
with rough leading edge (\’ 8mm, z 194 cm, 0 34 mm). 


In the conditions of the present experiments relating to the boundary 
layer, with transition by preturbulence, we have attempted a similar 
comparison, with a view to a possible extension of Taylor’s hypothesis. 

Figure 21 presents, for \’ 6 = 0-24, the comparison of the two auto- 
correlation curves R,,(7, X, = Ys = X3 = 0) (at positions s = 74cm and 
> = 79cm) with the longitudinal correlation curve R,,(7 = 0,.X,,.V. = 3 = 0) 
obtained along a line .Y, = 0, which here is not very different from the 
streamline. In this case there is still a satisfactory concordance between 
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the autocorrelation and the longitudinal space correlation, established up to 
X,/M ot about 2-5 (i.e. X,/V = 5:2 ms), but for X,/M > 2 (or X,/V >4ms), 
the correlations being smaller than 0-1, the comparison is no longer 
conclusive. 
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Figure 19, Transverse space-time correlations in the boundary layer on a flat plate 
with rough leading edge (v’ = 24mm, =z = 194m, 6 = 34mm). 
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Table 2. 

Further verification, including the measurement of the longitudinal and 
transverse space-time correlations R,,(7, X,, X, = 0, X3) along the line of 
maximum correlation, with an optimum delay T = 77, (see below), shows 
(figure 30) that this correlation retains high values of the order of 0-5 up 
to the distances (Y,/5), — 4; or (X,/M)p— »; given in table 2, but these 
distances decrease rapidly with the distance y’/6 of the upstream wire from 
the wall. 
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Measurements of the longitudinal and transverse space-time correlations 
Ri (7, X,, = 0, have been made in the case of the transition by 
preturbulence, using the old type of probes. ‘The upstream hot wire A 
was fixed, at a distance z = 74 cm from the leading edge, and at a distance 
y’ from the wall. The hot wire B was placed at the distance 1, downstream 
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Figure 20. Optimum delays 7 — T; for transverse space-time correlations in the 
boundary layer on a flat plate with rough leading edge. 
of A, Its distance y to the wall and consequently X3, was varied. The 
delay 7 was also varied. The results are presented in figures 22 to 25, 

which correspond to the following conditions : 


figure 22, y'/8=0-06, =1mm, X, = 25-4 mm; 
figure 23, vd = 0-06, y’ = 1 mm, X, = 50°83 mm; 
figure 24, y’ [5 = 0-24, y’ =4+mn, X, = 25-4 mm; 
figure 25, = 0-24, = 4mn, X, = 50-8 mm. 
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‘They follow the results published (Favre, Gaviglio & Dumas 1954b) for 
the conditions = 0-71, v’ = 13 mm, Y,=49mm._ The correlation 
is found to reach a maximum for an optimum value 7), of the delay, which 
is a function of .Y;. 
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Figure 21. Check on Taylor’s hypothesis for the boundary layer on a flat plate 
downstream of a grid — 4mm, = 0-24). 


Figure 26 gives the values of 7',, and their comparison with the computed 
values. ‘The computation is made as previously (Favre, Gaviglio & Dumas 
1954b), that is, by assuming (as a first approximation in the case of these 
experiments) that the optimum delays T,,, differ from the delays T, by time 
intervals which are equal to the delays T,. compensating for the transport due 
to the mean flow* ; 

T,, ~ T,;+T,. 
This relation is still satisfactory when, the upstream wire being at a distance 
from the wall y’/6 = 0-24, the downstream wire itself is at a distance 
v6 > 0-04; also when y’/d = 0-06 for y/5 > 0-09. For shorter distances, 
the relation is less satisfactory, either on account of the disturbances brought 
about by the probe or because supplementary terms should be introduced. 

Figures 27 and 28 indicate that, as we have shown elsewhere (Favre, 
Gaviglio & Dumas 1954b), the introduction of the optimum delay 7°, 
increases the correlation in a large proportion (of the order of 2:55 to 5-5 
in the cases examined). ‘The maximum maximorum of the space-time 
correlations, Rmax(7' = T,,,, X2 = 0, X3), corresponding to given is 
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* is given, to a sufficient approximation, by = 2X, (V,+-V>,). 
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Figure 23. Transverse and longitudinal space-time correlations in the boundary 
layer on a flat plate downstream of a grid (yv’ = 1mm, y’/d = 0-06, 
X, = 508mm). I. First series of measurements. II. Second series of 
measurements. 
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Figure 24. Transverse and longitudinal space-time correlations in the boundary 
layer on a flat plate downstream of a grid (\" 
X, = 25-4 mm). 
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Figure 25. Transverse and longitudinal space-time correlations in the boundary 
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Figure 26. Optimum delays T = 7,,, for transverse and longitudinal correlations 
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Figure 27. Comparison between transverse and longitudinal correlations with zero 
delay (J = 0) and transverse and longitudinal correlations with optimum 
delay (JT = T,,). Transition by preturbulence (z = 74 cm). 
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obtained for a certain value of Y,. Likewise the maximum of the space-time 
correlations, Rymax(7T = 0, Y,, X, = 0, X3), corresponding to given Xj, is 
obtained for 7’ = 0 and for a certain value of X;, ditferent from the previous 
me. In figure 29 have been drawn the lines of maximum space-time 


R 
| 
4 
4 
6 ~4 2 0 2 4 X 
3 
Figure 28. Comparison between transverse and longitudinal space-time correlations 
with zero delay (7 — 0) and transverse and longitudinal space-time correlations 
with optimum delays (7 7). Transition by preturbulence (z = 74 cm.). 


correlations with optimum delay (7 = 7,,), and also the lines of maximum 
correlations with zero delay (7 = 0). Outside the boundary layer, these 
lines coincide with the mean streamlines (Favre, Gaviglio & Dumas 1953 b), 
but inside the boundary layer, they are markedly different from the 
streamlines. ‘lable 3 summarizes the results. 
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Figure 30 shows the space-time correlations outside the boundary layer 
(y'/6 > 1) and inside the boundary layer (y'/6 = 0-71, 0-24, 0-06), along 
the line of maximum space-time correlations, relative to the upstream wire 
position. S,, and S, represent the curvilinear abscissae of the upstream 
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| | (62 
| 50°% 2.2 | 0-080 
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0-71 49-0 | 0-160 
| 
| | | 


Table 3. 


y /S 
Figure 29. Lines of maximum correlation with optimum delay (J = T,,) —O— 
Lines of maximum correlation with zero delay (JT — 0) -—-—- x) -———- 
Stream-lines (downstream of a grid). 


and downstream wires along this line. When 7'= 7., one gets the curves 
of space-time correlation along the line of maximum space-time correlation, 
with optimum delay. As we have noted above these curves show that the 
correlation retains high values of the order of 0-5 up to distances (.X1/5)p — 9.5 
which decrease rapidly with the distance y’/6 of the upstream wire from 
the wall. 
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Space-time correlations in a turbulent boundary layer 


In figure 31 are reproduced the space-time isocorrelation curves with 


optimum delay 7 = T,,, for the distances to the wall of the upstream wire 


y'/6 = 0-06, 0-24 and 0-71. 


Measurements are being continued, with improved accuracy. 
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Figure 31. Space-time isocorrelation with optimum delay, in the boundary laver 


on a flat downstream of a grid. 

‘These researches were made at the Laboratoire de Mécanique de 
l’Atmosphére (I.M.F. Marseille) for the Office National d’Etudes et de 
Recherches Aéronautiques (O.N.E.R.A.) with the aid of the Ministére de 
lAir and of the Centre National de la Recherche Scientifique. 
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On the hydrodynamic stability of two viscous incom- 
pressible fluids in parallel uniform shearing motion 
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SUMMARY 

A new problem in hydrodynamic stability is investigated. 
Given two contiguous viscous incompressible fluids, the fluid 
on one side of the plane interface being bounded by a solid wall 
and that on the other side being unbounded, the problem is to 
determine the hydrodynamic stability when the fluids are in 
steady unidirectional motion, parallel to the interface, with 
uniform rate of shear in each fluid. ‘The mathematical analysis, 
based on small disturbance theory, leads to a characteristic value 
problem in a system of two linear ordinary differential equations. 
The essential dimensionless parameters that appear in the present 
problem are the viscosity ratio m, the density ratio r, the Froude 
number F, and the Weber number W,, as well as the parameters 2, 
R (which is proportional here to the flow rate of the inner fluid) 
and c, that occur in the study of hydrodynamic stability of a single 
fluid. ‘The results obtained are presented graphically for most 
fluid combinations of possible interest. ‘The neutral stability 
curve in the (x, R)-plane is single-looped, as in the boundary layer 
case. ‘The calculated critical Reynolds numbers are higher than 
the values observed in liquid film cooling experiments. (In these 
experiments, the outer fluid is usually a turbulent gas, in which 
the thickness of the laminar sublayer is of the same order of 
magnitude as the liquid film thickness.) General agreement 
between the theoretical and experimental values exists for all 
critical quantities except the Reynolds number. Gravity and 
surface tension are found here to have a destabilizing effect on the 
flow, in agreement with experimental evidence. Semi-infinite 
plane Couette flow is a special case of the present problem and 
the known stability of this flow is recovered. The linear velocity 
profile of two adjacent fluids with the same viscosity, but different 
densities, is shown to be unstable for high enough Reynolds 
numbers. ‘The Reynolds stress distribution for a neutral 
oscillation in the general case is discussed qualitatively. 


I. INTRODUCTION 


The interest in the present problem stems from the desire to understand 
what happens when a liquid film flows over a flat surface, dragged along 
by a high speed gas. For certain liquid flow rates, the film surface becomes 
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wavy, and detached parcels of liquid from the main body are entrained by 
the gas and carried downstream. ‘This situation arises, among numerous 
engineering applications, in connection with film cooling of a solid boundary, 
as reported by Knuth (1954). 

A classical problem related to the one mentioned above is the generation 
of ocean waves by wind. ‘This problem has recently been investigated 
theoretically by Lock (1954), who was the first to have included in the 
analysis all the physical properties of the air and water. His calculations 
are incomplete, and the results obtained are quite different from the usual 
single-fluid results. 

The problem of the stability,of stratified motion of ditferent fluids has 
been studied by ‘Vaylor (1931) and Goldstein (1931), who did not include 
viscosity in their analyses. ‘Taylor investigated continuous and dis- 
continuous density and velocity distributions. Goldstein treated similar 
problems, his investigations being a generalization to heterogeneous 
stratified fluids of Ravleigh’s (1887) work on the homogeneous case. 
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Figure 1. Undisturbed velocity profile to be investigated. 


‘The problem of liquid film stability has been investigated experimentally 
by Kinney, Abramson & Sloop (1952) and by Knuth (1954), who were 
concerned with liquid film cooling applications where the gas stream was 
always turbulent. York, Stubbs & Teck (1953) studied the mechanism 
of disintegration of liquid sheets experimentally, and they proposed an 
inviscid model based on an extension of Lamb’s work (1932). 

The aim of this investigation is to solve the hydrodynamic stability 
problem when two viscous incompressible fluids, in two-dimensional, 
laminar, uniform shearing motion (figure 1) are perturbed by a small 
arbitrary disturbance. One of the fluids, from now on called the liquid, 
is bounded in the direction normal to the flow by a solid wall and by the 
second fluid, called the gas, of semi-infinite extent in the direction normal 
to the flow. 
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By ‘solving the problem’ is meant finding the relationships satisfied 
by certain parameters when neutral stability exists. ‘This implies that a 
neutral-stability hypersurface, whose coordinates are the physical variables 
involved, could be constructed that would separate regions of stability and 
instability. In other words, if the physical properties of the fluids are 
given, we wish to determine the minimum critical Reynolds number* at 
which instability begins. 

It will be realized from the above remarks that the chosen model is an 
idealization of the physical situation. ‘The most serious criticism is that 
the gas motion, in the actual case, is turbulent in most cases of interest, 
and the velocity profile is not a linear function of the distance away from 
the wall. Some justification for the present approach lies in the fact that, 
at least in the laminar sublayer of the turbulent flow, the gas is laminar 
with an almost linear velocity profile. (‘The ratio of gas laminar sublayer 
thickness to liquid film thickness in existing experiments is of the order 
of unity.) ‘There is still one other reason for the approach used (besides the 
obvious one of greater tractability when the turbulence is ignored), for 
it has been definitely shown by Zondek & ‘Thomas (1953) that semi-infinite 
plane Couette flow is always stable. Now, this situation is a special case 
of the present model, and it is interesting in this connection to know how 
a discontinuity in density or viscosity affects the stability of uniform shearing 


motion. 


II. ‘THE BOUNDARY VALUE PROBLEM AND ITS SOLUTION 


1. The Orr-Sommerfeld differential equation and its general solution for plane 
Couette flow 

‘The task now is to formulate mathematically the problem of stability 
of two-dimensional laminar motion. ‘This has been done in the past by 
numerous authors (e.g. Lin 1945). For the sake of completeness, however, 
a brief description of the derivation of the disturbance equation will be 
given here. 

Let all coordinates and velocities be made dimensionless by the use 
of a reference length 5, and a reference velocity Ue Consider a basic flow 
in the x-direction with a velocity profile U(y). The Navier-Stokes and 
continuity equations in the xy-plane can be perturbed by assuming the 
velocities in the v- and y-directions and the pressure to have the form 

U(v) + u(x, v, t), z(x,y, t), P+ p(x, v, t), (1) 
where the lower-case symbols indicate small quantities, and ¢ is the 
dimensionless time (i.e. time x U,/5). The introduction of (1) into the 
Navier-Stokes equations leads to two linear partial differential equations 
from which the pressure can be eliminated by cross differentiation and 
subtraction. ‘The result is a linear partial differential equation containing 

* Defined, for given flow physical properties, as the smallest value of R for 
which incipient disturbances will become amplified, i.e. the smallest value of R of 
the neutral stability curve in the (x, R)-plane. 


F.M. 2A 


: 
ky 

4 

iss 


346 S. Feldman 


u and wv as the dependent variables, with x, y and ¢ as the independent 
variables. ‘he equation of continuity guarantees the existence of a stream 
function A(x, y,¢), such that 


ow ous 
==, =, 
CV OX 


and which, when used in the differential equation last mentioned, reduces 
it to a linear partial differential equation in terms of only one dependent 
variable y(x,y,t). The tact that this equation is linear is very important, 
since this means that superposition of solutions is allowable. ‘This implies 
that if any arbitrary disturbance is decomposed into its Fourier components, 
it is then sufficient to solve the problem for one general sinusoidal oscillation. 
After the solution is obtained, it will be necessary to consider all possible 
frequencies and see how they affect the behaviour of the solution. 
In order to separate the variables in the partial differential equation 
for yb, let 
Jt) = (2) 
where x is the wave number, assumed positive without any loss of generality, 
and ¢ is the complex wave velocity which may be expressed as 
C= C,+1¢;, (3) 
where ¢, is the wave velocity and c; allows for amplification of disturbances 
if ¢, > 0, damping of disturbances if c; < 0, and neutral disturbances 


if ¢ =0. ‘The partial differential equation for y% then reduces to the 


Orr Sommerteld equation 


where the primes indicate derivatives with respect to y, and KR is the 
Reynolds number pdU,/j: (p is the density of the fluid and y the viscosity). 

If the velocity profile of the undisturbed flow is a linear function of y, 
as in the case of interest, 0” = 0, and (4) becomes 


— 224) = — +249), 0) 


i 
which is a linear fourth-order total differential equation in the complex 
y-plane. As first pointed out by Lin (1945), equation (5) has four linearly 
independent solutions, which are analytic functions of y and entire functions 
of the parameters c, x and aR. We will now solve (5). A transformation 
which was first pointed out by Orr (1906) and later independently by 
Sommerfeld (1909) ts 


b'(y)— = (6) 
which when used in (5) yields 
—[iaR(U—c) = 0. (7) 


Now, let 


so that (7) can be written as 


h’(z)+ 2h(z) = 0, (9) 
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which is the so-called Stokes equation. ‘The solution can be obtained in 
terms of contour integrals by using Laplace’s method (Morse & Feshbach 
1953, pp. 582-585). The result is 

h(z) = 2kyh,(z) + 2kyh,(2), (10) 
where the factor 2 is introduced for convenience in later calculations, 
k., and k, are arbitrary constants, and 


h(x) = exp(st + dt, hk = (12)"*e-™*, (11) 
hz) = exp(2t + 423) dt, = (12)!%e— (12) 


in which L, and L, are the contours of integration shown in figure 2, 
h, and h, are entire functions of z, and ¢ is a complex variable of integration. 
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Figure 2. Paths for contour integrals in the solution of Stokes’ equation. 


‘The functions ,(z) and h,(z) may also be written in terms of Hankel functions 
of order one-third as 

It will be recalled (see, for instance, Copson (1935)) that the Hankel functions 
are of oscillatory nature, H{!}, the function of the first kind, being damped 
exponentially as |3232| becomes large, while H{?} increases exponentially 
under the same conditions. 

From (8) and (10) we have 
C(y) = G(y) + 


where L(y) = &,(2), = h.(2), (14) 
which when inserted in (6) gives ; 
‘he complete solution of (5) is then 
H(¥) = bi(y) + Ro boly) + + ba(y), (16) 
2A2 


} 
; 
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where = o{y) = e", 


1 


bf) | dt — 07] M(t) dt. 


_ 


tis a variable of integration and should not be confused with the notation 
used to denote time. ¢, and ¢, are given by (14). 


2. Differential equations for the present problem 
When the phenomenon under investigation involves the  stratitied 
motion of two fluids, an equation like (5) must be used for each fluid, with 
the proper boundary conditions at the interface. Let 6 and U, as used to 
render (5) dimensionless be, respectively, the height of the liquid sheet 
and its surface velocity. Lower-case and capital Greek letters will represent 
conditions in the liquid and gas respectively, and the subscripts / and vg 
will denote quantities evaluated in the liquid and in the gas. We then 
choose 
as our representation of the disturbance stream function in the liquid 
and gas respectively. ‘lhe disturbance velocities then become 
ay 


ex 


(19) 


while the corresponding Orr Sommerfeld equations are 


(U,— ed" 234) = — — + 24), (<y<1), (20) 


— 220) = — — 202" +04), (1 <y < w), (21) 
( = 
8U, 
where R, = U, = ¥. (22) 


Equating the shear stresses at the interface of the basic flow gives 


U,=1+ (23) 
Mg 
‘The solutions of (20) and (21) can now be written, from (16), as 
= Ry bi(¥) + Ro holy) + + bly), (O<y<1), (24) 


(vy) = K,®,(y)+ K,®,(y)+ < &), (25) 


where the ¢’s are given by (17), and the ®’s can also be obtained from the 
same equation by using the appropriate value of = in the gas, given by (8). 
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Boundary conditions 

The boundary conditions will be written in terms of dimensionless 
variables. 

Let subscripts 1, 2 and 3 used with the coordinate y denote, respectively, 
the wall, interface and infinity, i.e. y, = 0, vy. = 1 and ys; = x. 

At the wall (y = vy, = 0), both components of the disturbance velocities 
must vanish; 1.e. 

= 9, = 0. (26) 
At the interface (y = y. = 1), the following conditions must hold: 
(a) Both fluids move together with no vacuum layer between them; 


thus 
v V2) vy(¥2) = 0. 27) 
(6) ‘There is no slip between the fluids in the direction of flow; thus 
— = 0. (28) 
(c) ‘The shear stress must be continuous; thus 


(d) Because of surface tension, the normal stress is discontinuous 
according to the relation 
= —Fyyr 
where p, is the ettective pressure caused by surface tension. If o@ is the 
surface tension coefficient and /, the radius of curvature of the interface, 
p, can be written as 


= 


where = (3) | 


the negative sign having been chosen because p, is to be positive when 
d*y dx? is negative. ‘The normal force aren can be rewritten as 


Py — pit + py = 9. (30) 
At infinity (y = v3 = 0%), the dinvarbente must vanish; i.e. 
Vg) = 0. (31) 


Equations (26) to (31) are the eight boundary conditions necessary for 
solving the system of two-fourth order differential equations given by (24) 
and (25). 


+. The secular equation 

Equations (13), (14) and (17) show that when v-> x, ®,(y) and 
D(v)-- «. In order to satisfy the boundary conditions given in (31), 
it is necessary that 1n (25) 


K, = K,=9, (32) 


1 

1 

a 

) 

| 

or 
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which when used with (19) allows the boundary conditions (26) to (31) to 
be rewritten as 


ky $y, + ko ho, + + hy by = 9, (33) 
hy bi + Robo + bis + hada: = 0 (34) 
ky + Ro boo + hg + Ry — Ky — Ky = 0, (35) 
Ry biz + Re boo + + Ry — K, = 0, (36) 


+ + hyo) + + + 


m = 


4 


a? 
a? 


where $4, = $23 ete., and 
R=R; F= U3/gs; W=p,dU3/c. (39) 
R, does not enter in the relationships of this section because it is clear, 
from (25) and (26), that 

R, = 


For the details of the derivation of (38), reference may be made to the 
author’s thesis (Feldman 1955, Appendix <A). 

Thus, the essential dimensionless parameters of the problem have 
been defined as x, R, c, m,r, F and W, where F and W denote, respectively, 
the Froude and Weber numbers. It is seen that Rk, has disappeared from 
the problem as an explicit parameter. ‘This seems to be in line with 
experimental results (Knuth 1954, pp. 362-363), which show that the 
inception point of unstable disturbances is independent of the gas-stream 
Reynolds number. 

If the preceding set of equations, (33) to (38), is to have a non-trivial 
solution for the k’s and K’s, the following relation, the so-called secular 


- m+ —-— |Rrm®,, —|[(1 —c)Rrm —i3a]®,, — - =0, (33) 
m F(1—c) 
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equation, must hold: 


by, bo 0 0 

by 0 0 

hie fre fe tre — Fy m — 

S12 §22 832 S42 — 
where 
faz = n= 1,2,3,4; 
n= 1,3; 

E + * 1] ~ ¢)R n= 1,2,3,4; 


= E + | arm. 


—[(1—c)Rrm — — w= 


5. Solutions of the boundary value problem 
‘The purpose of our calculation from now on will be to find the function 
that the secular equation represents. 
From (22), (23) and (38), it is clear that in the liquid 
=. U’'=1, (<y<1}),} 
and in the gas } (+1) 
U=1+m(y-1); U'=m 0). J 
‘The six solutions involved are, from (8), (11), (14), (17), (24), (25) and (32), 
$(v) =e", (42) 
=e", (43) 


¢,(v) = | —t,)\€(t,) dt, — | exp| — (45) 
1 J 


| 
L 


a 


/ 
dt, — | exp} — dty>, (44) 


®(y) (46) 
®,(y) = | dt, — | — — ty) 
(+7) 
s(t) = = = | exp(z)t + dt, (48) 
C,(t)) = = — | exp(2,t + 47°) dt, (49) 
k - 
= = | expleyt + 40%) dt, (30) 
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w 
bo 


= \ 1(xRrm*) 1) |}. (52) 


and for clarity, instead of using y in the liquid and gas, the variables ¢, and ¢, 
were introduced. In (44) and (45), the lower limit of integration was taken 
as y = | for convenience in later calculations. (For an alternative refer to 
discussion following (55).) Later on, many expressions will have to be 
evaluated at y= y, = 1. With the choice made here, the calculation is 
simplified since several integrals vanish, i.e. $52, $52, Pag and 
dys. In(47) the lower limit was taken as infinity, because two of the boundary 
conditions required the solutions to vanish there: no other choice would 
have been satisfactory. ‘he functions needed in (40) will now be written down : 


= 1, =i, 
| 
dy» Pog = Dio 
“0 ) 
bs exp( — xt))C,(t;) dt; — | exp(at)¢,(t) dt, 
~ 1 


= | exp( — xt))C,(t,) dt, | exp(xt))C,(t,) dt, | 
J 1 


d,. = 0, day = 0, 


deo = 2C,(t; = 1) = | 
bi, = = 1) = 
(53) 
re “0 
| exp( — dt, + | exp(xt,)Co(t,) dt, 


| 1 
{ 


dig = 204i, = 1) = 
Pag = = 1) = 


| | x(t, —1)}¢,(t,) dt, | exp{x(t, — 1)}E,(t,) dt,, ! 
| exp{ —a(t,—1)}¢,(¢,) dt, + | exp( x(t, — 1) dt 
= 20,(t, = 1) + = 2h,[2,(1)] + 


= 26i(t, = 1) + = 3, 
Here 31) and z,(1) denote the values for t, = 1 and ¢t, = 1 respectively. 


| 
| 
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We shall next change the determinant of (40) into a form more 

convenient for calculation. The procedure is as follows: 

(a) Divide the last row by R. 

(b) Multiply the 3rd, 4th and 6th columns by («R)!?/'¢5., (aR)!?/ dys 
and (xRrm?)'?/':, respectively, where the accent to the left of a 
symbol means the leading term in the asymptotic expansion (aR > 1) 
of the corresponding function. 

(c) From (60), insert all the zeros for the terms that vanish. 

(d) Multiply the 5th column by —1, and rearrange the columns so that 
column 5 becomes 3, 3 becomes 4, and + becomes 5. 

The result is then 


| A, 0 b, 

| Ay V, 0 ay b, 0 

tA, Vs Fy 6, 

| V, 0) 0 

| . 


where the first three columns involve the inviscid terms (except for the 
last row, which has terms of O(1/R) and will be neglected in comparison 
with the terms kept in the calculation), and the last three columns involve 
the viscous solutions. ‘lhe term‘ viscous solutions’ refers to any of ¢3, $4, 
‘),, or their derivatives, regardless of where they are evaluated. ‘The 
meaning of tne svmbols in (54) is given by the following relationships: 


(1— 
3=%2, %Vy= doe, V5 = 


a, = “b32, = (35) 
3 
by = = (AR) | (59) 
hb, = (aR)!  b, = | 


L 

} 

4 

3 
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x 


(60) 


/'Doo 
szf 


It is worth, remarking that the lower limit for the integrals in (44) and (45) 
could alternatively have been taken as zero instead of unity. ‘This would 
have lead to a fourth-order determinant instead of the sixth-order one 
in (54). Each element of the alternative determinant would have been not 
only more complicated than the one of (54), but the clear distinction between 
viscid and inviscid solutions (of which use is made for the calculations of 
the neutral stability lines) would not have been possible. 


Behaviour of the viscous solutions for aR > | 

‘Thus far, no restrictions of any kind have been made in the analysis. 
We will restrict the discussion to the case xR > 1, this being sufficient 
to solve the problem. One of the major tasks is the evaluation of the 
integrals in (53). ‘The only previous discussion of a similar integral is 
the one by Hopt (1914), who assumed «Rk < 1; thus his discus jon 1s 
not suited to the present study. In order to find out how th. .ecutral 
stability curve (i.e. for given m, r, F and W) behaves for large wala s of 
the parameter «R, it will be sufficient to keep only the terms «| highest 
order in «R in the asymptotic expansion of the pertinent functinrs. As 
will be shown later (following (71)), the case of interest for neutral stability 
(ct. (3)), is the one where ¢ = c, — 1. A detailed description of the method 
used for obtaining the integrals, the order of magnitude of the errors 
involved, and the calculation of ¢,, as an example of the procedure used, 
is given in Appendix B of the author’s thesis. If xR > 1, from (51) and (52) 
it may be assumed that 
= 2(0) = 3c ~ —i(aR)!3c, 
3), = 3(1) = —c) ~ i(aR)' 31-0), 


Zo = 234 (1 —c)m ~ —c)m, 
(61) 
Now, let 
e+ = exp[3(aR)! 2e©44(1 — c)3?], | 


exp[ — 3(aR)! — c)32], > (62) 
E = — — | 
so that (58), (59) and (60) become 
a, = a,,(aR) Of(2R)'4e ], 
= + Of[(aR)' te], 
= Ag, + O(a2R) 
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= by + + O[(aR)Me>], ] 
b, = b,, + aR) + O[(aR)! 4e+], [ 
¢ 
5 = *+ bs» + O[(aR) | (64) 
= boy + + + | 
65 
where 
| 
a, % 34] = —e (66) 
! 
a;, = 1, Ay = | 
b.. — _ Sinha 
| 
a9) Cosh x 3/2 
— gil = —e-“ ~c) | 
| 
be, = e 7A) be (67) 
bo = a» %Sinhz, boo = 
22 
— C 
m \32 m \32 | 
l—c l—c 
/ 12 
| 
/ m 2 (rm? 1/2 
Cy = tele “(rm*) C1 = Fa 
| 
| 


The secular equation (54) can be expanded in terths of a sum of products 
of the viscid and inviscid terms. Equations (62) to (68) are then used for 
the purpose of determining the important terms to be kept in the final 
result. 


i 
€ | 
) 
4 
on 
Bers 
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7. The eigenvalue problem for zR > x 
It can be shown that for xR > «, the secular equation reduces to 
F(a, c,r,m, F, W)=0, 


(1 —c)a[ g(a, 7, m) +r] — +1- mr) = 0), (69) 


or 


which, when solved for c, gives* 
(a? W+1- rm) + W+1—rm)* F} 


9(a,1r,m) +r] 


(70) 


there pt. 12 
g(a, m) = | 1] (71) 


sinhz  tanhz m/|r)' 


the radicand being always positive, i.e. c is always real. ‘This means that 
for R-- x, the flow ts neutrally stablet. Also, since for cases of interest 
rm —_ 1, c ts always less than unity as R-- «, ‘This is quite an important 
conclusion, since for finite Reynolds numbers, the computation is different, 
depending on whether c is less or more than unity. It will therefore be 
assumed, in all the work for finite Reynolds number, that c < 1. The 
final calculations will bear out this assumption. 

A different approach to the case of infinite Reynolds number would 
be to neglect viscosity at the outset, 1.e. in the Orr- Sommerfeld equations. 
The differential equations then become of second order, and by relaxing the 
proper boundary conditions, the problem could again be solved. This 
has been done, and the equation obtained for ¢ is the same as (70) except 
that the function g(%,7,m) is not given by (71) but by 


2(a,r,m) = coths, (72) 


which agrees with (71) only when mr > 1. 

The limiting processes used to obtain (71) and (72) were ditterent. 
‘Theretore it is not surprising that these results are not identical. From 
a physical point ot view, the disagreement is inacceptable. Although no 
physical explanation of the discrepancy has been found as yet, the dis- 
agreement 1s immaterial trom the present standpoint because for the fluid 
of interest, viz. a liquid gas combination, m/r > 1 and the results agree. 


S. The case of finite Reynolds number 
‘The fundamental equations for the neutral stability curve in the xR-plane, 
for a given physical situation (i.c. for fixed gas-liquid density ratio r, 


“The root with the negative sign in cont of the radical can be proven to be an 
extraneous root introduced when solving (69). 

+ It is worth remarking that for the case of ‘ inviscid Couette flow ’ of a single 
fluid the Orr-Sommerfeld equation has no non-trivial solution, while for two layers 
of fluids there 1s only this one non-trivial solution for the eigenvalue problem. 


. 


| 
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liquid--gas viscosity ratio m, Froude number F, and Weber number W), are 
H (x, 
F(a, c)’ 
each point on the curve having a particular wave velocity. ‘Che functions 
F (a,c), Ga, c) and #%(x,c) are defined by 

F (4,0) = x + [e (Ag — V4) +27, sinh 


r 


G(a,c)— (a,c) = 0, (aR)? = — 


(74) 


2a 5 17 1 


m\ 12 1 
[e (A, — V,)+27, sinh a] E + 41 (=) ] + 


[((V,—Agle +27, cosh — 


2a = 5 m {m\'2sinhz 
~ | = (3) cosh + 


1 sinh « 
2a(1 (l-c)!* « (9) 


H (4, ¢) = —292(1 —c) + af(A,+ 27, cosh 


— 2x( —A,)cosh (76) 


a? 1 
1, =| m+ fre (79) 


The solution of (73) for c as a function of « has to be obtained numerically*, 
because of the impossibility of solving explicitly for c or « in terms of the 
other. Once a pair of values of x and c is known, straightforward calculation 
leads to the value of («R)'?, it being possible to evaluate R immediately. 

A digression is permissible at this point in order to point out the meaning 
ot Reynolds number R~ R, in the present problem. Since the velocity 
profile is a linear function of distance, for given liquid physical properties 
R is proportional to the liquid flow rate, and therefore is a constant once 
the flow rate is specified. 

* All the numerical work was carried out on an IBM card-programmed electronic 
calculator using an 8-digit floating decimal system. 
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Now that a method for determining the neutral stability curve has been 
described, it is necessary to decide which region, on either side of it, is 
stable or unstable. ‘The calculations indicate, as is reasonable to expect, 
that for a disturbance of a given wavelength, the flow is stable for Reynolds 
numbers smaller than a Reynolds number &,,, corresponding to a neutral 
disturbance; and similarly the flow is unstable for R > R 
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Figure 3. Neutral stability curves: wave number vs liquid Reynolds number. 
Gravity and surface tension forces are neglected. Liquid—gas viscosity 
ratio 10. 


The results obtained from the present analysis, of which figure 3 is a 
typical example with F = W = 0, indicate that the shape of the neutral 
stability curves in the «R-plane is similar to the boundary-layer case*. 
The effect of varying the gas-liquid density ratio r could be stabilizing or 
destabilizing, while increasing the liquid-gas viscosity ratio m always 
tends to stabilize the flow, provided m is not small compared with unity. 


*'The reason for obtaining only one solution to the problem, and not two as 
Lock (1954) found, is that in the present case each fluid is stable by itself, and the 
only reason for the existence of instability is the discontinuity in the physical properties 
of the fluids. In Lock’s case each fluid is probably also unstable. On the other 
hand, in the present treatment instability regions which do not extend to R= & 
have not necessarily been eliminated ; this may occur when c > 1. The analysis 
would have to be extended to this case in order to investigate this possibility. 
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For «0, the equation for the lower branch of the neutral stability 


curve 1s 
xR = const., (80) 


as can be derived analytically. ‘The upper branch seems to have a nearly 
horizontal tangent for R > oo; this was obtained numerically. 

The influences of non-vanishing gravity* or surface tension forces on 
the above results should be alike, since they enter as a sum in the secular 
relation (40). Figure 4 shows that for large Froude F and Weber numbers WW, 
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Figure 4. Effect of gravity and surface tension forces on neutral stability curves. 


gravity and surface tension destabilize the flow. Obviously, surface tension 
effects are negligible when the disturbance frequency is small, but become 
very important at large frequencies. ‘The effect on the neutral stability line 
is in this case to raise the upper branch so that its asymptote for x > 1 is 

a/R = const., (81) 
which can also be obtained analytically. 


* The gravitational field was assumed to act downwards in the flow configuration 
indicated in figure 1. 
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Figure 6. Rate of change of amplification factor with Reynolds number as a function 
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All the effects discussed thus far concern the stability of the flow, and 
the direction (left or right) in which the neutral stability curve moves in 
the aR-plane when the physical quantities are varied. 

There are two or more items to be discussed before the results of our 
numerical calculations are summarized. ‘They are (a) the behaviour of the 
magnitude of the wave velocity; (6) the amplification or damping rate of 
disturbances in the neighbourhood of the neutral stability line. 

With respect to (a), it is enough to mention that all waves travel at 
speeds less than the velocity of the liquid-gas intertace, a typical case 
being given in figure 5. 

On the neutral stability curve, the imaginary part c; of the complex 
velocity ¢ vanishes. It is possible, nevertheless, to compute the rate of 
change of c; with respect to Reynolds number. ‘lhe important result is 
that this is always a positive quantity which, as a function of « (figure 6), 
has a peak near the critical Reynolds number. ‘This means that the curves 
¢; = const., in the «R-plane, would be packed close together when in the 
neighbourhood of the critical values of x and R. 

Since the quantities of greatest interest are the parameters corresponding 
to the critical value of the Reynolds number, the large number of calculations 
which we have made may suitably be summarized as follows. 


9. Values of critical quantities and discussion of results. Comparison with 
experiments 

Since the physical case of interest is the one where gravity and 
surface tension forces are small, detailed calculations were carried out 
for F=W= <x. ‘The critical values are presented in figures 7 to 10, 
from which the following facts can be gathered: 

(1) As m and r — 1, the flow is completely stable. 

(2) For a given gas-liquid density ratio +, an increase in the liquid gas 
viscosity ratio m always increases the stability, and decreases the amplification 
or damping in the region away from the neutral stability curve. 

(3) For very small r, the flow is completely stable. 

(4) For a given m, there is always a value of r for which the flow is most 
unstable. 


(5) As r — 0, the wave number « tends to a value of 0-6 approximately, 
and the wave velocity c becomes 0-1 approximately. 

(6) For the special case of air and water at a temperature of 100°C 
and a pressure of one atmosphere (m = 10, r = 0001), the value of the 
critical Reynolds number is 60000. 

(7) The critical conditions depend for given physical properties only 
on the liquid flow rate. 
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The significance of the results quoted as items (1) and (2) can be 
understood in terms of the known universal stability of plane Couette flow 
(from now on abbreviated as P.C.F.) between walls of arbitrary spacing. 
An explanation follows. 

In the special case when the two fluids have the same density and 
viscosity (equivalent to the case of a single fluid), the velocity profile becomes 
a single straight line. Item (1) shows the flow then to be universally stable, 
in agreement with the known result for P.C.F.; a check on the analysis 
is thus obtained. 
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Figure 7. Critical Reynolds number. 


When the result quoted as item (2) is interpreted for the limiting casc 
ot a very viscous liquid, the motion is again always stable. The very viscous 
liquid could just as well be considered as a solid and this case again reduces 
to PALF, 

From item (2) and the above discussion, a result can be deduced which 
has not been obtained directly by any calculation. ‘This is that, for an 
arbitrary density ratio, the flow is completely stabilized when m -> 0. 
The reason for stabilization is that this limiting case of the flow occurs 
when the gas becomes so viscous that it could be replaced by a solid wall, 
which again reduces to P.C.F. between two walls at a finite spacing. 
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‘Therefore, in figure 7 there would exist a curve for some small value of m 
that would be farthest to the left, and for smaller values of m the curves 
would again be displaced more and more to the right as the liquid-gas 
viscosity ratio is decreased. 
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Figure 8. Critical wave number. 


Items (5) and (6) will now be compared with some experimental 
observations. Most experiments with liquid films have been carried out 
in horizontal round tubes where the liquid flows along the inner surface 
dragged by a high-speed turbulent gas. It should be kept in mind that 
while the gas layer is of infinite thickness in the theoretical model, the 
ratio of laminar sublayer thickness in the gas to liquid film thickness in the 
experiments was of the order of unity. 

The experimental neutral wavelength A observed by Knuth (1954) -for 
all liquid flow rates was about 10 film thicknesses. Considering that 
g = 27/A, Knuth’s findings check with the critical value of item (5). Since 
the wavelengths in the neighbourhood of the critical value are the most 
amplified, these might be the only ones visible in an experiment. Knuth 
might have observed these values, which, although seemingly neutral, 
could have been slightly amplified. 
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Figure 10. Critical values of the rate of change of 


S. Feldman 


= t 
m=} 
T fj 
= 
T + 
| | 
r / 
10 | 
= 
T 
| | 
| 
| I > | 
| | 
— 
| / 
= | 
I/ | 
-3 | 
10 + ‘3 
| 
| 
| | 
| | } j 
F=W=00 
-4 | 
ile] 
10? 10" 
Corit 
Figure 9. Critical wave velocity. 
| 
rh WA 
| | | 
| 


ac, 


number. 


amplification factor with Reynolds 


j 
; | 
| 
10 
} | 
| 
r 
10 
| | 
| 
10°} 
| 
-4 
10 3 <8 <4 
10 10 10 10 10 


Stability of two viscous fluids in uniform shearing motion 365 


Item (6) implies that for a liquid film 0-005 in. thick, the critical 
liquid-gas interface velocity is 100 ft/sec, which is one order of magnitude 
larger than the values obtained from experiments in liquid films. This 
discrepancy would seem to indicate that the observed instability is not 
simply laminar instability of uniform shearing motion. ‘The fact that the 
computed critical Reynolds number is much higher than the experimental 
values could be due to the fact that the velocity prefiles of both fluids in 
the analysis were assumed to be straight lines. If curved profiles were used 
it is conceivable that the critical Reynolds number could decrease. ‘This 
possibility is suggested by the change in the value of critical Reynolds 
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Figure 11. Influence of gravity and surface tension forces on critical Reynolds 
number. 


number when going from plane Couette (R= *”) to plane Poiseuille 
How (R.,;, = 11560, based on the maximum velocity and width of the 
channel). Curvature in the velocity profile of the liquid could exist in 
the case of laminar flow when there is a pressure gradient in the flow 
direction. One of the things that has not been accounted for in the present 
or suggested analysis is the effect of turbulence of the gas stream, which 
also possibly influences the stability of the flow. 
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Gravity and surface tension have a destabilizing effect (see figures 11 
and 12), as remarked in the discussion in §8. ‘This influence of surface 
tension was observed by Kinney, Abramson & Sloop (1952, p. 10), the 
relative change in the critical Reynolds number being of the same order 
of magnitude (within a factor of at most 2) as the value found here 
analytically. 

‘ Item (7) agrees with the experimental fact reported by Knuth that the 
inception point of instability is a function of the liquid flow rate, and 
independent of the gas flow rate. 
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Figure 12. Influence of gravity and surface tension forces on critical rate of change 
of amplification factor with Reynolds number. 


Before concluding, it will be helpful to try to gain some physical insight 
into the stability problem by looking at the energy of the disturbed motion. 
This will be done in Part III. 


10. Accuracy of calculations 


The accuracy of our calculations, based upon the viscous functions 
given in §6, Part II, is restricted by the fact that xR was assumed large. 
This was interpreted as meaning that in (61) the right-hand members are 
good approximations to the left-hand members, and permitted the use of 


‘ 
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asymptotic methods. ‘Terms of O[(aR)~'] were neglected when compared 
with terms of O[(aR)~!*] in the secular equation. A check on the assumption 
mentioned was done for each calculation of a point on the neutral stability 
curve. The result of this checking showed that the assumption was incorrect 
for small values of r, shown dotted on figure 7 only. Nevertheless, the 
trends in that region are probably correct. ‘This is inferred by comparing 
the dotted sections with the curve for m = 50, which is valid over almost 
all of the regions shown: they seem to form a reasonable famuty. 

For small values of r, the last of equations (61) for z,, should be 
approximated by the first term and not by the second, as for larger values 
of r. If this innovation were made, a valid calculation for the dotted part 
of figure 7 would then be possible. 


11]. ‘THE REYNOLDS SHEARING STRESS 
A different way of looking at the stability problem, due to Lin (1954), 
consists of following the history of the disturbance energy, which, for 
damping or amplification, changes as a result of the action of the Reynolds 
shearing stress 7,*. It is then enlightening to know its distribution across 
the stream. 
Foote & Lin (1950) have shown that 


— exp}2x¢; — ¢¢’), (82) 
and d¢'), (83) 


where wz is the distance or time average of the product uz, and ¢ and ¢’ 
dicate complex conjugates of ¢ and ¢’. 
At the interface, for the viscous case, 


= (34): 


(Tal Pog = (85) 
The inviscid case 

The amplitude functions ¢ and © for the inviscid case are real functions. 

Therefore (82) shows that 

— = const. = 0, (86) 
from which it follows that the Reynolds stresses are zero across the flow. 
This means that there is no mechanism for transferring energy between 
the basic flow and the disturbance, i.e. any disturbance will just subsist, 
without damping or amplification. 

As can be seen, the method of this section is extremely useful, since 
the important result of §7 regarding stability has now been rederived in 
a few lines without any calculation. 

* The work done per unit volume per unit time by the basic flow is tdU/dy 


(Schlichting 1950), and converts energy from the basic flow into the disturbance 
when dU/dy > 0. 
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The viscous case 

We will start the study of this case by showing that the Reynolds stress 
is continuous across the layer where U(y) =c,. The Orr-Sommerfeld 
equation (4) can be rewritten as (¢ is complex) 


— 22h” + (87) 


U" i 


which, on multiplying by 4, subtracting its complex conjugate and regrouping 
vields 


U—c2 2R U—c 


_ bb — + 


dy 


(88) 


I;quation (88) can now be introduced into (83), and since UL” = 0, we have 


dy 4R 
which can be integrated across the layer where l= c, and remembering 
that in our case U = y, 


(ar. c+0 (7), c—0 = 


lim - 4 _ [log « — log( — + 


Since loge —log(—e) = +1n(2n+1), 


equation (90), with ¢, = 0, gives for the jump in Reynolds stress [7,] across 
the laver where U = c 


[7] = + bb — + ][ + (YD) 


where a particular » should be chosen for the branch of the logarithm 
being used. (Note that the quantity in the first bracket on the right-hand 
side of (91) is real, and the quantity in the second bracket is purely imaginary.) 
Since [7,] must be real, and p +R + 0, we have at v = ¢ 

bd + — + 6d) + = 0. (92) 


Lin (1954) has shown that the Reynolds stress in a viscous fluid grows 
positively and very rapidly with distance away from the wall in a very thin 
laver, and then stays about constant. 

We now have enough information to build the complete picture of stress 
distribution for the viscous case. Starting at the wall, the Reynolds stress 
is zero, and as we proceed outwards it grows at a rapid rate. It then levels 
otf and since viscosity is not very important, (86) will be almost satisfied 
(i.e. tr, = const.). There is no jump across the layer where U = c. As we 
get to the liquid gas interface, the stress is discontinuous, the ratio of the 
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values on both sides of the discontinuity being given by (85). ‘The stress 
on the gas side cannot be zero, since by (84) this would mean that the stress 
in the liquid is zero. Far away in the gas stream, the Reynolds stress must 
vanish. Figure 13 shows qualitatively this stress distribution. 

Lin (1954) found an important relationship in the theory of hydrodynamic 
stability in a very simple way. ‘The principle used by him was to equate 
the stress in the fluid adjacent to the wall computed by two different methods, 
starting from the wall in one case, and from the main stream in the other. 
‘This method has not yet been made successful here, the reason being that 
the magnitude of the jump of stress across the interface is unknown, only 
the ratio being known. In order to compute the jump, one would have 
first to compute the stress at the interface, on the gas side, by calculating 
the function ®. ‘This would destroy the simplicity of Lin’s method. 
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Figure 13. Distribution of Reynolds stress for a neutral oscillation. 


Therefore, looking at the Reynolds stress in the present case was not 
as fruitful as in Lin’s case. It would nevertheless be interesting to find the 
quantitative distribution of stress for a self-excited disturbance and for a 
neutral disturbance. ‘This would show whether most of the energy input 
into the disturbance comes from the gas or the liquid. From the functions 
presented in this paper, it would be possible to calculate the Reynolds 
stress for a neutral disturbance. ‘The case of a self-excited disturbance 
is much more complicated, and the necessary amount of numerical work, 
as envisaged at the present time, is prohibitive. 


IV. CONCLUDING REMARKS 
The discrepancy between experiments on the stability of liquid films 
and the present theory is confined to the fact that the theoretical value of 
the critical Reynolds number is larger than the experimental one, all other 
quantities being in agreement. ‘The results obtained here seem to indicate 
that a study of the hydrodynamic stability of two fluids with curved velocity 
profiles would be worth-while. Such a study might definitely settle the 
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question whether or not the large-scale disturbances observed in liquid-film 
cooling experiments are due to laminar instability. 

The distribution of the Reynolds stress across the stream remains a 
problem of great importance to the physical understanding of hydrodynamic 
stability in the present example. 

' ‘The model chosen for the analysis, although the simplest possible, has 
yielded a number of new and interesting results. ‘The most important 
conclusion to be drawn from this investigation is that a discontinuity of 
viscosity or density has a destabilizing etfect on uniform shearing motion. 


This paper is based on the thesis submitted by the author in partial 
fulfilment of the requirements for the degree of Doctor of Philosophy at 
the California Institute of ‘Technology, June 1955. ‘The author wishes to 
express his gratitude to Professor W. D. Rannie, for suggesting this problem 
and for his critical guidance during the course of this investigation, and to 
Professor C. C. Lin for his very helpful suggestions and personal interest. 
The author is also grateful to the Daniel and Florence Guggenheim 
Foundation and to the California Institute of ‘Technology for making 
possible his continuation of graduate study through the grants, respectively, 
of a Fellowship in Jet Propulsion and a Graduate Tuition Scholarship 
during the 1953-54 and 1954-55 academic vears. 
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Oscillations of tide and surge in an estuary of 
finite length 


By J. PROUDMAN 


Edgemoor, Verwood, Dorset 
(Received 18 February 1957) 


SUMMARY 

This paper concerns a narrow basin of uniform cross-section 
open to the sea at one end and closed at the other. An incident long 
wave of prescribed general form is supposed to enter from the sea 
and to represent the combination of tide and surge as generated in 
the sea. ‘The solution of the linear terms of the equations of 
continuity and motion gives the reflection of this wave at the head 
of the estuary. ‘This paper gives the next approximation when the 
non-linear terms are retained, the second-order motion being 
made determinate by the condition that, at the mouth, it reduces 
to a progression towards the sea. 

The chief results relate to the surface elevation at the head of 
the estuary. When the first order elevation there increases steadily 
to a maximum, the effect of the ‘shallow water terms’ is to make 
high water higher and earlier, while the effect of the ‘frictional 
term’ is to make high water lower and later. For a short estuary, 
the interaction of the tide on a surge, due to a given sequence of 
meteorological conditions over the sea, is to make it higher when 
its maximum occurs at the time of tidal high water than when its 
maximum occurs at the time of tidal low water. This is directly 
opposite to the corresponding result when the estuary is of 
infinite length. 

INTRODUCTION 

In two papers (Proudman 1955 a, b) I have discussed the dynamics of a 
progressive wave of tide and surge in an estuary, and particularly the 
interaction between the tide and the surge. For equal sequences of 
meteorological conditions over the sea, | showed that the apparent height 
of a surge whose maximum occurs near to the time of tidal high water 
is less than that of a surge at the same place whose maximum occurs near 
to the time of tidal low water. But the tides of the ‘Thames Estuary, for 
example, constitute a standing oscillation much more nearly than they 
constitute a progressive wave. 

In the first of the papers, I also gave an approximate solution of the 
differential equations relating to an estuary with ‘a barrier, but without 
using end-conditions which would completely determine the motion. In 
1956 Doodson considered the problem presented by a gulf, and gave 
numerical solutions for the cases of a tide prescribed at the mouth and of 
a surge prescribed at the head. Because of reflection from the head, it is 
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inappropriate to prescribe a surge at the mouth of an estuary. Also, any 
attempt to deal with discontinuities of basin at the mouth, involves 
reflections from the mouth, both on the seaward side and on the estuarial 
side. 

In this paper, a prescribed incident wave is taken to progress up the 
estuary, and the remaining part of the motion is taken to reduce, at the 
mouth, to a wave travelling down the estuary. ‘The incident wave is taken 
to represent the combined tide and surge which are generated in the sea, 
and the remainder of the motion is taken to represent the origin of the 
wave which, in the open sea, will diverge away from the mouth of the 
estuary. 

In order to avoid a greater mathematical complication, | suppose the 
cross-section of the estuary to be uniform. I follow the same mathematical 
method as in my earlier paper, except that I now use the differential 
equation tor the current, whereas previously I used that for the elevation 
of the water-surface. Because of the more precise end-conditions of the 
present paper, the results are more definite than those of that part of the 
earlier paper which relates to an estuary with a barrier. 

It appears that the surface-elevation at the head of the estuary depends 
on the primary elevat.on there at the same time and for a previous interval 
during which a progressive wave could travel twice the length of the estuary. 
When, during this interval, the primary wave increases steadily to a 
maximum, the tendency of the shallow water terms of the differential 
equations is to make high water higher and earlier, while the tendency of 
the frictional term is to make high water lower and later. For a short 
estuary, the current in the estuary, and the time and height of high water 
at its head, are approximately independent of friction; while the interaction 
of the tide on a surge, due to a given sequence of meteorological conditions 
over the sea, is to make it higher when its maximum occurs at the time of 
tidal high water than when its maximum occurs at the time of tidal low 
water. ‘This result is directly opposite to the corresponding result for a 


progressive wave. 


2. NOTATION AND DIFFERENTIAL EQUATIONS 
Denote by : 
¢ the acceleration of gravity, 
h the undisturbed depth of the water, supposed uniform, 
a the length of the estuary, 
v distance down the estuary, x = 0 being at the head and v = a at 
the mouth, 
t the time, 
€ the elevation of the water-surface, 
u the current down the estuary, 
k a numerical coefiicrent of friction, which will be taken as 0-0025 ; 


and write c= (phy. 
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Take also x 
E=t+-, =, (1) 
Cc 
so that 
t= 23+), <== (2) 
ct cn CX c€ cn 


Accents will be used to denote derivatives of functions with respect 
to their arguments. Suffixes 1 and 2 will be used respectively to denote 
terms which are of the first and second orders in the ratio of the primary 
surface-elevation to the depth of water. 

The equation of continuity is 


cl 


- 
= (h+ Qu} + = 0, (5) 
CX ct 
and the equation of motion is 
Cu Cu are k (6) 
cl Ox ° dx h 


while uw = 0 where v = 0. 
‘The product terms on the left-hand sides of (5), (6) will be called the 
‘shallow water terms’ of these equations. 


3. GENERAL FORM OF SOLUTION 


‘he elimination of ¢ from the first-order terms of (5) and (6) gives 


crt? ox® 2 exct h ct 


and the general solution of the first-order part of this equation, which 
makes u, = 0 where x = 0, is 


: x x 
—F(t+ 2) (8) 


where F() denotes any function which is physically interpretable. It 
follows from the first-order part of (5) that 


so that » ( 
(t+ (9) 


no term independent of ¢ being required. 
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‘The terms of (8), (9) in F(t+./c) are taken to represent the prescribed 
incident wave. The solution of the fundamental equations, as far as the 
second order, will be of the form 


+F(t-—-}+, 10) 
h ( h ( 
—-= (11) 
c 
The terms of (10), (11) in F(t—./c) represent a wave travelling down the 
estuary, and, to satisfy the condition prescribed in $1, so must the terms 
in at x =a. The condition for this is 
Us 


where v = a, and this will be used as the determining condition in this paper. 


4+. SECOND ORDER TERMS 
For the second approximation, the equations (5), (7) lead to 


h—+— = 13 
Ay a9 A9 A 
O"Us O"Us 1 u-) k 
of © Ax? >1 i) Ox at h ot ( ) 


respectively, and the transformation (3), (4) applied to (14) gives 


Now (9), (8) may be written respectively as 


>1 


= F()+ F(a), (16) 
“t= — F(g)+ F(a), (17) 


and substitution from these into (15) gives 


$ 


‘The general solution of (18) is of the form 


(19) 


where 


4+—— = Gy — 
( 
c& en h 
c& cn cn h cy 
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and X, ¥ are arbitrary functions of €, 7 respectively. One integration of 
each of these equations gives 


k 
+ FOF) + FO) Fa), 


kc 
Pa) FE) + FE) (FO) — FO]: 


no additional arbitrary functions being required. Further integrations 
give 


$= FAO! + FO | " 


FO) [FO FO] (20) 


hj; 
b= — | F(8) db + 
ke 


the additive functions of € and 7 respectively in these integrals being s¢ 
chosen that ¢ = ys = 0 where €—7 = 0. 

From (19), (20), (21), and on taking 


A= —4f(), 4f(m), 


it follows that 


= + - + FD} | a0 - 


1 ke 
FO- FO) FO) FO- Fa) Fey} 
+l). (22) 


In (22), us = 0 where €—y = 0, whatever the function f( ) may be. 
Transformation of (22) back to the independent variables x, ?¢, by using 
the equations (1), (2), gives 


| -|F®)- F(t F(t =) |} dé — 
=) + f(1-*); (23) 


12 
| 
Ly 
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and then substitution from (8), (9), (23) into (13) gives 


~ 
—~_ 
~ 
ale 
| 
ale 
+ 
| 
all 
t 
—~ 
~ 
| 
ale 
+ 


f(t *) *), (24) 


since the terms due to differentiation with respect to the upper and lower 
bounds of the integrals in (23) balance. Integration of (24) gives 


Substitution from 


(2 
\ \2 


3), (25) into the condition (1 


for all values ot ¢, so that 


—~ 
te 
~~ 
~ 
bol 
~ 
~" 
| 
| 
— 
D> 
a 
| 


f(t) = - F(t) - P(e - 


| Cc J 
Lhe | ( x) | 
x 
(25 
( c} ( (£9) 
| 
a\ {3a = a _ 
(1 5 | F(#) — 
c ¢ a ac J 
1 ke a | ( a 
*) ral *) | d 
| 
a F(t)— [F(t)— F(@)] dé. (26) 
+h J t-2a¢ 
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Substitution from (26) into (25) giv 


) gi 
= 


t 


Fle Ble bol vo 


t—x/e—2a 
+t—a2/e 
| Fl] F(t *)- Fe | dé + 
t+2/e—2a/c 
1] ke rt+2/c x 


and a similar formula may be obtained for w/c. 
For the second approximation to be a valid approximation, it is 
necessary for |¢, to be small compared with the maximum value of ¢. 


5. ELEVATION AT HEAD OF ESTUARY 
At the head of the estuary, where x = 0, it follows from (10), (27) that 


; = —}F(t)F(t—2a/c) + 


+t 
Cc 2 

F(t) — F(t) — F(@)] dé. (28) 

It is then seen that the elevation of the water-surface at the head of the 

estuary depends upon the first-order elevation there at the same time and 

for a previous interval during which a first-order progressive wave could 
travel twice the length of the estuary. 

When the first-order elevation at the head of the estuary has a positive 

maximum at ¢ = 0, then F’(0) = 0 and F'’(0) is negative. ‘The formula (28) 

then gives, for that time, 


2F(0)+ 2(0) 2a/c)- | F(0)F(—2ac)— 
= F(0) — F(@) [F(0)— F(@)] dé, (29) 


= —!F(-2a F(-—2a/c)} + 
+ Re — F(—2aJe) {F(0)— F(—2aJc)}, (30) 
F.M. 2C 


x 2a 1 x 2a 

‘ 

1 x x 1 x 2a 

1 x x 2a 3a-x “\.. x 

4 c c a 

+f — F(t--)F(t--)- F(t+-) | F(0) dé — 

x *t-+2z/c 

—§F(t--) | F(@) dé 

4 

AB 
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and (30) is of the second order in F. 


Now ec 7 (=) (=) 
et et 1-0 Ot? 


to the first order in ¢, and, as high water occurs when c¢ ct = 0, it follows 
from (30) that the time of high water is given, to the first order in F, by 


F(—2ac)} + | F(#) dé — 


F(0) — F(—2a/c) | F(O)— F(—2a/c)}. (31) 


Also, cl 
C= (¢) 


so that, to the second order in F, the height of high water is given by (29). 
Suppose now that 


0 < F(—2a‘c) < < F(0), for —2a/e <6 <0, 


“ F(0) 


bol 


so that F’(—2a_c) is positive. Then from (29), (31) it can be shown that 
the tendency of the shallow water terms of the fundamental equations is 
to make high water higher and earlier, and that the tendency of the frictional 
term is to make high water lower and later. 
When there is no surge and the first-order tide is harmonic of amplitude 
2hA and period 27 « so that 
F(t) = Acosot, 


then the formulae (29), (31) for the height and time of high water at the 
head of the estuary become respectively 


1 2ca_ 1 toa 
= 24+ A*-— — -cos - ~cos — — 
. 20a 1. 
— ~— — — 2sin— + -sin— }]}, 
Af 30a 1. 40a lAkef, 4 20a 4ca\ 
In the tidal example taken by Doodson, A = 128 ft, 27/¢ = 12 hours, 


aajc = 27/5, 2hA = 5 ft, so that A = 0-0195, ke/ho = 6-9, } Akce/ho = 0-034. 
These figures indicate that the method of the present paper should give 
an approximation to the motion. But for an estuary, the depths are smaller, 
so that the values of ke/ho are larger, and the values of A may be much 
larger. In these circumstances the approximation will only be valid if aa/c 
be small. 
6. SHORT ESTUARY 

Now suppose that, during the time taken by a progressive wave to travel 
the length of the estuary, the elevation at any one place changes by only a 
small fraction of its maximum value. Then aF’/c will be small compared 
with the maximum value of F. 


= 
5 
: 
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On expanding the second-order terms in ascending powers of x/c as 
far as x3/c#, it follows that 
| dé = | (F(t) +(6—t) F(t) + F'(t)} dé 
t—wie ~t—Z/e 


x ix 


and similarly it may be shown that 


rt 
| = F(t) F(t) + F’(t). 


/ t-2ac 


Also 


= F(t) F(t)| dé 


~ 


and similarly it may be shown that 
| F(8)— F(t— xc) #)— dd = F(t), 
t—2z/e 
F(8) [F(t)— F(8)] dé = 3 a 
It then follows, from (25), (23), (26) respectively, and on omitting the 
ss ania t of F and f, that 


3 
13 1x3 
+5 ( FF (33) 

f= — + 3 + FF. (4) 

From (10), (32), (33), (34), it then follows that 
(FF" +3F’F’)- 
k 


(35) 


lI 


+ FF )+ S(FF +4F'F"), (36) 


the argument of F being ¢t, except where otherwise shown. 
2C2 


Li 

_ 

F() F(d, 

: 

3 

ok 


conditions over the sea may be regarded as giving rise to equal functions S(t). 
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It will be seen that, to the order x? c®, the current is independent of 
friction. 
When F(O) > 0, F'(0) = 0, F’'(O) = 0, 


the formulae (29), (31), for the approximate height and time of high water 
at the head of the estuary, become 


> = 2F(0} 1— — F"(0)+ = F’(0)}, (37) 
h J 
2a a” F(O) { 2 a® 
t — F(0)+ F"(0)— 5 F*(0)}-— 5 = F'(0). (38 
Both these results are independent of friction. 
When again F(t) = Acosot 


A, o being constants, the tormulae (35), (36) become 


? 4 2 
h 2.4 cos — cos ot + A* | in Zot + (a* — x") 


( 
(a*+ 3x") cos 2ot + = —;(a’— x’) sinot sinot >, 
u 
24 sin — sinot + 
c 


? 5 
 L0°/ : t4o7ax 
4°. — + |sin Zot + —— COs Jct +, 


respectively. ‘The formulae (37), (38), for the height and time of high 
water at the head of the estuary, become 


| 
we! bo 
“ 
tw 
te 


respe ely. 

‘The validity of the expansions in x ¢ requires that cac < 1. For the 
above formulae to give a valid approximation to a solution of the fundamental 
equations, it is necessary for (ca ¢)A, (ko*a® he*)A to be small. 

As an example, take h 40 ft and a semi-diurnal tide. ‘Then 
¢ = 24:5 miles per hour, o = 27 12 hours™!, so that co = 47 miles. With 
a = 20 miles and 4 = }, this gives 

aa ko?a? 


~A=014, 


A = 0:39. 
c he? 


7. SEPARATION OF TIDE AND SURGE 
Take F(t) = T(t)+ S(t), (39) 


where 7() denotes a tide and S()asurge. Equal sequences of meteorological 


Since the timing of the meteorological conditions which generate surges 


a 
4 
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is independent of the timing of the astronomical forces which generate the 
tides, the frequency-distribution of the first-order surges S will be 
independent of the phases of the first order tides 7. 

Substitution from (39) into (35) with x = 0, gives for the height of the 
combination of tide and surge at the head of the estuary 


4 2a” 
T’+S'(T'+S’'), (40) 
the argument of 7’ and S being ¢. 
The predicted tide would be given by 
4a 2a* 4 ka’ 


so that the apparent surge is see ‘a subtracting (41) from (40). ‘The 
formula for the — surge is thus 


4 4 
= = 28+" Ss" +5 + 4 
2 4 
€ 
+ ka’ 


When the first-order surge has either a maximum or a minimum so 
that S’ = 0, the formula (42) reduces to 


2@ 2a a* 
é 


+3 


c 
772 
c? 
and this is free from frictional influence. 
When the terms in a*c? are negligible, and, at the same time, the 
first-order tide is either at high water or low water, so that 7” = 0, the 
formula (43), for the apparent surge at the head of the estuary, reduces to 


2 a2 
= 5") +TS’). (44) 


Now suppose that 
S > 0, 


so that the first-order surge has a positive maximum. When the first-order 
tide is at high water, 


T = A, T’ = -o°A, 


i 
43 
= 
? of 
i 
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hA denoting the amplitude of the tide and 27/c its period, and the apparent 
surge (44) is given by 


2 
= 28(1- 55") (45) 
When the first-order tide is at low water, 
T = -—A, = oA, 
and the apparent surge (44) is given by 
2 
j 28(1- 58") - (46) 


‘The excess of (45) over (46) is 


4a? 


and this is positive. 

It is thus seen that, at the head of a very short estuary, the ettect of the 
interaction of the tide on the surge is to make the apparent surge higher 
when its maximum occurs at the time of tidal high water than when its 
maximum occurs at the time of tidal low water. ‘This result, which is due 
to the shallow water terms of the differential equations, is directly opposite 
to the corresponding result for a progressive wave. ‘Vhat result, in the case 
of an estuary of uniform cross-section, is due to the frictional term. 

‘To examine the order of magnitude of the difference (47), take the 
particular case in which, when the first-order surge is at its maximum, 
S” = —o?S. ‘The difference in the heights of the apparent surges is then 
equal to the height of the first-order surge multiplied by (So7a?,c?)4. With 
the figures at the end of §6, this factor is 0-47. 
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SUMMARY 

An analogue method for solving certain quasi-linear hyperbolic 
partial differential equations is presented. ‘The characteristic 
lines are formed by scanning electronically the screen of a cathode- 
ray tube. ‘The boundary conditions are introduced in the form of 
an opaque mask. ‘The solution appears as a picture on the screen 
of a second cathode-ray tube. ‘he experiments demonstrate the 
feasibility of the approach, but the development of the machine 
has not been carried beyond this stage. 


GENERAL DESCRIPTION OF THE METHOD 

‘There are well-established analogue methods for solving Laplace’s 
equation, such as the use of an electrolytic tank or of a soap film, to mention 
just two. An extensive bibliography of such methods has been given by 
Higgins (1956). On the other hand, there has not been comparable progress 
in finding analogue methods for solving the wave equation. Step-by-step 
numerical or graphical solutions are available. For example, for two- 
dimensional flows that are hyperbolic everywhere, the method of 
characteristics is applicable. It is essential to recapitulate certain features 
of the method of characteristics for the present point of view. 

In the case of steady two-dimensional supersonic flow, there are two 
tamilies of characteristic lines which can be regarded as a curvilinear 
coordinate system. In the case of a uniform supersonic flow, these 
characteristic lines become two sets of parallel straight lines each forming 
the Mach angle with the flow direction. If the departures from uniform 
flow are small, these characteristics will differ only little from the straight 
‘undisturbed’ characteristics. 

For the linearized theory of supersonic flow, the characteristics are 
assumed to have fixed direction (independent of the solution), and only 
the higher approximations involve the ‘ warping’ of the characteristic lines. 

Scanning methods are essential features of the television technique, and 
they can be quite readily adapted to generate the characteristic lines. 
A spot moving on the screen of a cathode-ray tube can sweep according to 
a scanning programme of arbitrary complexity. Present commercial 
television uses a scan pattern that consists of one family of parallel straight 
lines (usually in the horizontal direction from left to right) which are formed 
sequentially in time. 
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If a scanning programme is devised so that the ordinary television scan 
is performed twice, irst tilted with a positive then with a negative angle, 
two families of fixed characteristics can be generated. For practical reasons 
it is more convenient to generate the two sets of lines as part of a single 
scan pattern as given in figure 1. By using symmetrical triangular waves 
for vertical deflection and saw-tooth waves for horizontal deflection with 
only slightly differing frequencies, one obtains a dense scan pattern. ‘The 
scan lines also appear to be ‘reflected’ from the top and bottom edges of 
the scanned area. 


VERTICAL 


HORIZONTAL 


Figure 1. Scan pattern. 


By choosing the light intensity of the cathode-ray screen as the dependent 
variable, solutions to the linear wave equation can be obtained by simply 
holding the light intensity constant during each scan stroke, but varying 
it in an arbitrary manner between scan strokes. 


‘The assignment of light intensity for each scan stroke depends on the 
boundary conditions, and a simple method can be devised for introducing 
such boundary conditions. Suppose the light intensity is the integral of a 
function which is zero in the field, but not on the boundary. ‘The spot will 


y y 
2525 | 

| 
| 
| 
82525 
t 


‘a 


An electronic analogue for supersonic flow 38 


have constant intensity after it leaves the boundary and this constant can 
be set at the boundary. Pictures were obtained with the use of two 
cathode-ray tubes as indicated in figure 2. Identical scan patterns were 
applied to the two tubes. ‘The first one, labelled SCANNER, sweeps in 
front of a cut-out mask giving light only when crossing the diamond-shaped 
slit. ‘The light intensity is picked up by a phototube, and is integrated by 
an integrating circuit. ‘The integrand is different from zero only when the 
scan spot passes the cut-out area. ‘lhe value of the integral is displayed on 
the second tube labelled DJSPLAY. ‘The light intensity is constant along 
scan lines outside of the ‘body’. Figure 3 (plate 1) shows the results of 
such an experiment. 


VLLL LLL DW 


SCANNER DISPLAY 


Figure 2. ‘'wo cathode-ray tubes. 


‘lhe integrator can either be reset to zero at the extreme values of the 
vertical scan (figure 3(a), plate 1), or permitted to hold its value and be 
reset only when returning to the left vertical edge (figure 3(b), plate 1). 
This gives a different set of ‘wall’ boundary conditions, 

In an actual flow problem, the boundary conditions are imposed by the 
slope of the surface and not by the integral of body thickness as in the above 
example. However, it is quite possible to construct masks so that the 
integrated values of light transmittance inside the body contour are 
proportional to the slopes of the contour. 

The analogue for the linearized supersonic flow problem as described 
above appears to be quite trivial, and the interest in this method was mainly 
due to the fact that a certain amount of non-linear behaviour is rather 
easily introduced. If the disturbances are small. but not infinitesimal, 
the principal effect of the non-linearity is to warp the characteristics by 
changing their local angle (and position) due to the change of the local 
Mach number and flow direction (see Lighthill 1954). In the analogue the 
simplest way to alter the characteristic angle is by changing one of the 
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scan velocities. For example, the horizontal velocity of the cathode-ray 
spot can be altered by an amount proportional to the magnitude of the 
disturbance. Since the scan is obtained by integration, this change in 
velocity also modifies the position of the scan line. 

For the non-linear operation, the boundary conditions are introduced 
in the same manner as for the linearized problem, namely, by scanning a 
mask as shown in figure 2. The video signal obtained from the phototube 


MASK SCANNED 


INTEGRAL 


HORIZONTAL SWEEP 


Figure +. Formation of the horizontal sweep for non-linear operation. 


by scanning the cut-out mask is integrated, and this integrated signal 1s 
considered as the ‘disturbance’ function. ‘This signal can be used to 
modify either or both the horizontal scan velocity and the light intensity of 
the cathode-ray spot on the DISPLAY tube. Figure + shows the modifi- 
cations of the horizontal scan by such an arrangement. Figure 5 shows 
a three-tone mask for a diamond aerofoil giving the correct boundary 
conditions when used with an integrator. Since the slope can be both 
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positive and negative, the background is a neutral grey, and there are both 
positive (light) and negative (dark) areas contributing to the integral that 
represents the slope of the surface. 


BODY 


4 
Z 


Figure 5. Proposed three-tone mask for introducing proper boundary conditions in 


case of diamond aerofoil. 


ANALYSIS 


1. The second-order equation 
‘The steady two-dimensional flow of a perfect gas has been the subject 
of extensive theoretical studies during the last hundred years. If the flow 
is assumed to be isentropic and the viscosity and heat conduction are 
neglected, the assumption of potential flow is justified. 
[et us introduce the following quantities : 
v, vy, Cartesian coordinates, 
d(x, v), velocity potential, 
u(x, y), 7(x,y), velocity components in the xe and y-directions, 
q(x, absolute magnitude of the velocity, 
, local speed of sound, 
g,, maximum possible velocity of the gas, 
y, ratio of the specific heats. 
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‘The local speed of sound depends only on the absolute magnitude of 
the local velocity, 


(1) 

12 
also om (1a) 
The governing equation of the flow is (see, for example, Howarth 1953) 
(c?— —26,$, br, = 0. (2) 


‘This second-order non-linear partial ditferential equation is designated as 
‘quasi-linear’ because the highest-order derivatives (¢,,,¢,,,,,) all 
occur linearly, and only lower-order derivatives occur non-linearly. 
‘The equation is elliptic where the flow is subsonic, 1.e. 
and hyperbolic where the flow is supersonic, 1.e. 
9 
‘There are tew exact solutions of this equation. 
Let us introduce the concept of an undisturbed flow and a superimposed 


small perturbation 


d(x, v) = U,x+e,4(x, y); (3) 
then u= 
and v= 

where ¢ = when +" = 0. It is also convenient to use M, = U, ¢,, 
the Mach number of the undisturbed flow. From equation (1) we find 
that 


By substituting ¢ from (4) and ¢ trom (3) into (2), we obtain a still exact 
equation for the disturbance potential ys: 


My + 12M, + = 9. (5) 
Equation (5) can be simplified by eliminating certain higher-order 
terms. Keeping only second-order terms and introducing 8? = 7-1, we 
obtain 


Dividing through by 1 —(y—1).7, ¢,, and again retaining only second-order 
terms in the expansion, we obtain 


= My (13), (7) 


yu 
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where « = [(y+1)+8}(y—1)]/28,. Equation (7), which is exact up to the 
second order, is in a convenient form for our present purpose. ‘The right- 
hand side is proportional to vw cv,cx, and it vanishes as a second-order 
quantity when the x-axis is turned into the direction of the local flow. 
Since the purpose of the present analysis is to demonstrate the effect of 
non-linearity, the right-hand side of (7) is neglected and a modified second- 
order equation used : 


= 9. (8) 

If we introduce 
e(x, Vv) = KM, (9) 
this becomes (8, = (10) 


2. The linearized equation 
The linearized equation for the present problem is the well-known 
wave equation 
Brads, %. (11) 
The general solution is 
(x,y) = F(wt+ B, vy) +G(x— v) (12) 
where F and G are arbitrary functions. 
The pressure disturbance is computed in the following manner. For 
polytropic changes we have 
ne 
pp * = const., (13) 
and the equation of state takes the form 
yp p = (14) 
We see that all thermodynamic properties are unique functions of the 
absolute velocity g. ‘The pressure takes the form 


= pl (15) 
In the linearized approximation, we have 

P= py — My) (16) 

py—p, Myth. (16a) 


If the wave equation is solved for a supersonic flow over a solid boundary 
y = f(x), the boundary conditions are imposed in the form that the flow 
direction must follow the tangent of the solid boundary; thus 
u dx M,+%,° M, 
It the oncoming flow is undisturbed, only one family of characteristics 
carries a signal (G = 0, or F = 0). One then obtains 

by = + Bits (18) 

(the sign depends on whether F or G vanishes), and 
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This shows that the local pressure disturbance is proportional to the local 
slope of the boundary. 

The lines €=x+8,y=const. and »=x-—f,y= const. are the 
characteristics, and they form a fixed network of oblique lines (Mach 
lines) along which the disturbances ‘ propagate’. ‘The linearized theory 
improves as the disturbance becomes smaller (very slender body), and is 
really valid only for the limiting case of infinitesimal disturbances. 

3. The analogue system 

‘Two families of lines are generated, and they are regarded as two sets 
of characteristics. With the use of a finite number of scan lines, there will 
be only discrete values of € and 7. 

The coordinates of the cathode-ray spot during the ath stroke are given 
as a function of time by 

X(t) = Xo +(t—t)X, (20) 

¥(t)=(t—ty)¥. (21) 
In general Y and ¥ are arbitrary functions of the time; Ay and ty are 
arbitrary constants. By making a total of \ strokes, the two sets of 
characteristics can be obtained in the following way. 

Choose for the nth stroke 


X,= €, form even, 
», for n odd, | 
t, = nll, (22 b) 
where At is the period of one scan stroke, also 
X= V(8,+«,) for even, 
V(8,+65,) for n odd, 
=(-1)'V. (22 d) 


Here ¢€, and 6, are numbers for 0 <n < \, and they are introduced by 


(22 c) 


the boundary conditions. For convenience, one can choose equal steps 
in €, and 7,, te. €&, = mp, and 7, =p. Using a dense set of lines, the 
subscript can be dropped, and £&,7 as well as 5,¢ can be considered as 
continuous variables. “hus 


(23a) 
y. 23 b) 


Observe that if « = 0 and 6 = 0, then € = const., and 7 = const. become 
the characteristics for the linearized problem. 
If e(x, ) and 6(x,y) are regarded as the dependent variables, it is easy 
to verify that they obey the first-order partial differential equations 
(B, +e)e, +, = 0, (24) 
(B, +8)5,—5, = 0. (25) 
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The variables e(x, y) and 6(x, vy) each also obey the following second-order 
equation : 

+5)*5,,—5,, = — 283 B, — 2883. (27) 
Naturally, these equations are valid separately, each for one family of 
simple waves only. By comparing (26) with (10), it is found that they agree 
for all except the last term on the right-hand side, but that is of the third 
order. This agreement can be interpreted in the following way. 

The analogue generates simple wave solutions of the second-order 
hyperbolic partial differential equations (26) or (27). They are also simple 
wave solutions of the simplified second-order equation for the flow problem 
if a third-order term is disregarded. ‘lhe dependent variable in both cases 
is the disturbance in the characteristic angle. For simple waves, the relation 
between e and the flow angle @ is monotonic. ‘This can be seen as follows. 

Using only the linearized approach from (18), we find (taking the upper 
half-plane) 


by By (28) 
The increment in the absolute magnitude of the velocity vector becomes 
dqiq = —d6/B,; (29) 
from (1), we get 

dec = —}(y—1)M? dq/q; (30) 

and, from the definition of Mach number, 
dM/M = dqq-dc c, (31) 
or dM/M = —(14+3(y—1)M?) d6/B,. (32) 


Also, the change in the cotangent of the characteristic angle with respect 
to flow direction becomes 


dB = — M*(1+ }(y— 1)M}?) d0/B;. (33) 

‘The characteristic makes an angle x with the x-axis given by 
cota = B, (34) 
= (35) 

After linearization, this becomes 
B, +e = B,+dB—M;? db, (36) 
= dB- do, (37) 

(38) 


Equation (38) indicates that the disturbance to be imposed at the beginning 
of each scan stroke must be proportional to the surface slope, and the 
ratio constant is a function of the nominal Mach number of the flow. 
Since the value of € is constant during a scan stroke, the simplest way 
of imposing its value at the boundary is by integration. An integrating 
circuit holds the value of « during the scan. It is reset after each stroke. 


392 Leshe S. G. Kovdsznay 


The signal is built up before the cathode-ray spot leaves the body contour 
The boundary condition can be introduced as a variable opacity or variable 
thickness of a zone within the body. 

Figure 5 shows one simple method for introducing the right boundary 


condition. 


4+. Light intensity of the screen 

The solution is obtained in the form of a photograph taken from the 
screen of the display cathode-ray tube. ‘The light intensity on the cathode- 
ray tube is however, not proportional to «(), but depends on both the 
luminosity of the spot and on the duration the spot spends in the local area. 
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Figure 6. Light intensity relations. 


In order to find the relations between € and the screen intensity, introduce 
H(t), the instantaneous luminous intensity of the cathode-ray spot, and 
I(x, v), the average luminous intensity per unit area of the screen. The 
intensity is held constant during a single stroke but varies from stroke to 
stroke. Note that 7(f) and /(x,v) do not have the same dimensions. 
Introduce for conventence 
i(t)7 (39) 


off 


ab 
where ah is the total scanned area and 7 is the time of the entire scan cycle. 
‘The local intensity obeys the equation of conservation of (luminous) energy 

I(x, y) Aw Ay = 1(t) > AT, (40) 
where ¥ \7 is the sum of time intervals the spot spends in the area Aw Ay. 
The luminous intensity /(v,\) can be expressed as shown in figure 6. The 
small quadrangle is bounded by €, = const., €,;+A€ = const., y, = const., 
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(a) 


Figure 3. Linear operation. (a) Integrator reset both at the end of each vertical and 
horizontal scan stroke. (6) Integrator reset only at the end of horizontal 
scan stroke. 


Figure 9. Non-linear analogue for supersonic flow. (a) Mask scanned. (4) Picture 
obtained on display. 
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Figure 8. Crowding of characteristics over concave surface. 


(a) 


Figure 10. Different degrees of non-linear effect. (a) Weak non-linear effect. 
(6) Strong non-linear effect. 
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and y,+Av=const. Let us assume successive scan lines correspond to 
equal increments in € (they are equidistant on the x-axis). The fraction 
of the number of scan lines that pass through the area is A€/2a (the factor 2 
is due to the fact that there are two families of scan lines during a full scan 
cycle). The fraction of the time the spot spends within the area during 
a single stroke is \y/d (the vertical velocity is constant), so that 


SAT _ Ay 
T 2ab 


Substituting this value in (40), we obtain 
21)(€) Ag = I(x, y) Ax, 


or = 


If the light intensity is modulated in intensity so that 


de 


the luminous intensity of the screen becomes 


This is a significant result since it was previously shown that € is analogous 
to pressure (or density) within the range of linear approximation. In 
wind-tunnel experiments, the xv derivative of the density is made visible 
optically by the schlieren method. ‘This is the most common visualization 
method used in routine wind-tunnel work. If the luminosity of the spot 
(intensity modulation, or ‘ Z-axis’) is varied according to the change of «€ 
from scan stroke to scan stroke, schlieren pictures result. 


EXPERIMENTS 


The experiments were carried out only to indicate feasibility of the 
method, and were not pursued to any degree of perfection. A flying-spot 
scanner previously described by Kovasznay & Joseph (1955) was modified 
for the present purpose, and many of the limitations of the original 
equipment determined the experiments that could be attempted. The 
block diagram of the equipment is shown in figure 7. 

The flying-spot generator produces a video signal corresponding to the 
transparent and opaque areas of the mask. A real image is formed in the 
plane of the mask by an objective lens and the signal obtained at the phototube 
is of the ‘yes or no’ type. A bi-stable trigger quantizes the phototube signal 
into a definite two-level signal (0 or 1). The reset integrator integrates the 
signal during one scan stroke. The integrator is reset to zero either at the 
end of the horizontal strokes or only at the end of both vertical and 
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and y,+Ay =const. Let us assume successive scan lines correspond to 
equal increments in € (they are equidistant on the x-axis). The fraction 
of the number of scan lines that pass through the area is A€/2a (the factor 2 
is due to the fact that there are two families of scan lines during a full scan 
cycle). The fraction of the time the spot spends within the area during 
a single stroke is Ay 4 (the vertical velocity is constant), so that 


AE Ay 
2ab 


Substituting this value in (40), we obtain 


or 


the luminous intensity of the screen becomes 


Kx, = (43) 


This is a significant result since it was previously shown that € is analogous 
to pressure (or density) within the range of linear approximation. In 
wind-tunnel experiments, the x derivative of the density is made visible 
optically by the schlieren method. ‘This is the most common visualization 
method used in routine wind-tunnel work. If the luminosity of the spot 
(intensity modulation, or ‘ Z-axis’) is varied according to the change of « 
trom scan stroke to scan stroke, schlieren pictures result. 


EXPERIMENTS 


‘The experiments were carried out only to indicate feasibility of the 
method, and were not pursued to any degree of perfection. A flying-spot 
scanner previously described by Kovasznay & Joseph (1955) was moditied 
for the present purpose, and many of the limitations of the original 
equipment determined the experiments that could be attempted. ‘The 
block diagram of the equipment is shown in figure 7. 

The flying-spot generator produces a video signal corresponding to the 
transparent and opaque areas of the mask. A real image is formed in the 
plane of the mask by an objective lens and the signal obtained at the phototube 
is of the ‘yes or no’ type. A bi-stable trigger quantizes the phototube signal 
into a definite two-level signal (0 or 1). ‘The reset integrator integrates the 
signal during one scan stroke. ‘The integrator is reset to zero either at the 
end of the horizontal strokes or only at the end of both vertical and 
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horizontal strokes. In the former case, the waves‘ reflect’ from the vertical 
edges of the scan pattern (see figure 3, plate 1); in the latter they do aot. 
If the masks were cut in such a way that the transparent portion along the 
scan line corresponded to the slope of the boundary, the signal obtained 
after integration becomes proportional to the ¢ disturbance. In the analogue 
for the linearized problem, we can simply modulate the light intensity 
with this signal. Figure 3 (plate 1) was obtained in such a way (with the 
switch, shown in figure 7, set on INTENSITY CONTROL). For non- 
linear operation, however, the horizontal scan velocity must be modified 
according to «. ‘This is obtained by the variable velocity sweep. This 
circuit is essentially a second reset integrator. In absence of a velocity 
control signal, it integrates a constant and X = 8, V, and it is reset at equal 
intervals. ‘The velocity control signal is added to the integrand, forming 
X = VB,+Ve as required by (22 c). 


FLYING SPOT SCANNER 


SQUARE WAVES | 
INTEGRATED PHOTOTUBE 


| 
[_ 1 CH 


MASK 
| . 

| VARIABLE say | 
— VELOCITY 
| | CAMERA | 
| | X (t) A] ! INTENSITY CONTROL | 
| | (FOR LINEAR OPERATION ONLY) | 

| | SWITCH | 
(NONLINEAR OPERATION ONLY) | 
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RESET PULSES 


== 


Figure 7. Block diagram of equipment. 


In the present experiments, the light intensity J, was not varied 
according to (42), but was kept constant. In this way the flow pictures do 
not actually correspond to schlieren pictures. 

A more serious defect of the present experimental set-up is that the: 
boundary conditions were not introduced correctly. ‘The masks were 
simply cut out of black paper, thus giving two levels of the video signal. 
Every smooth closed body contour, however, has both positive and negative 
slopes. In order to form both positive and negative disturbances, one needs 
at least a three-level discriminator instead of the bi-stable trigger, and this 
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was not available at the time of the experiments. ‘There is really one 
exception to this comment. Figure 8 (plate 2) represents a concave * body’, 
and the slope is approximately proportional to the local thickness so in this 
particular example the proper boundary conditions were applied. (A wedge 
with exponential contour satisfies the requirement that thickness is 
proportional to slope.) 

‘To modulate the light intensity in accordance with (42) would have 
required the use of a memory circuit to retain the value of € from the previous 
scan stroke, hence to form de/d& = e(€,,) —€(€,,_,). 

It should be emphasized that both of these improvements (proper 
boundary conditions by use of three-level discriminator and light-intensity 
modulation according to (42)) require only standard electronic design, and 
were not attempted only for lack of time available for the work. 

Figure 8 (plate 2) clearly shows the crowding of characteristics into the 
formation of an envelope. Figure 9 (plate 1) shows the photographic print 
of the cut-out mask and the ‘ flow picture’ thus obtained. Since the period 
of scan cycle is constant and the horizontal scan velocity is variable, the 
right edge of the picture is irregular as the total w distance is varying with 
the disturbance. In order to render the pictures more ‘life-like’, a small 
artifice was introduced. On the display cathode-ray tube intensity grid, 
a linear combination of the video signal and its time derivative was imposed. 
The negative value of the video signal from the phototube produced darkness 
in the interior of the body and the time derivative of the video signal gave 
a sharp outline of the body contours. Without this artifice the body contours 
would not even have been visible. 

Figure 10 (plate 2) shows two objects. In the picture in figure 10 (a) 
the non-linearity is ‘weak’, having low amplification of the disturbance 
signal. In the other picture, the non-linear effect is strong. 

As a possible extension of the method, the analogue for axisymmetric 
supersonic flows seems to be quite feasible. A more speculative proposition 
is the use of iteration. Using iterative processes, the first flow picture 
would be obtained by scanning a cut-out mask representing the boundary 
conditions. ‘The display would be photographed and would be introduced 
into the scanner to obtain a second approximation. In this way, the 
interaction between the two families of characteristics would be introduced. 

The extension to the one-dimensional unsteady problem, where the 
cathode-ray face would represent the x,¢ plane, is also obvious. It is 
interesting to note that in this case the analogue becomes exact for y = 3, 
as shown in the Appendix. 


The present work was carried out in the National Bureau of Standards, 
Washington, D.C., under a cooperative programme, of basic instrumental 
research and development while the author was a consultant. ‘The author 
acknowledges the enthusiastic help of Mr H. M. Joseph and Mr N. Newman. 
‘Thanks are due to Dr B. 'T. Chu of Brown University for contributing the 
special case given in the Appendix. 
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APPENDIX 
The one-dimensional unsteady flow of a compressible gas presented on 
an (x,t)-plane has a qualitative similarity to the steady two-dimensional 
flow. "The governing equations for an isentropic perfect gas are: 


cp cp eu 
Conservation,of mass, — +p= = 0. (A 1) 
at OX Ox 
Conservation of momentum, — +u— + = 0. (A2) 
pox 


By eliminating p, one obtains the pair of equations 


2 2 
t t (u t x. (« = 0, (A3) 


\ 


? ? 
> (« (u (« (A 4) 


Introduce the characteristic lines € = const. and 7 = const., i.e. with 
E = &(x,t), 7 = 7(x,t), so that 


dé = dx —(u+c) dt, dy = dx —(u—c) dt. (A5) 
Along these lines the characteristic relations are obtained: 
Along & = const., u+2c (y—1) = const. 
Along 7 = const., u—2c/(y—1) = const. 


The characteristic relations can be conveniently written as 
c+ = €(€), c—}(y—1)u = +4(y), (A 6) 
where e(€) and 6(y) are arbitrary functions of their arguments. In general 
the characteristics will be curved lines. However, if y = 3, 
u+c = ¢y+e(€), = (A7) 
Consequently, the equations (A5) can be independently integrated, giving, 
for instance, 
‘These equations represent two ‘ non-interacting’ sets of characteristics. 
Unfortunately all gases have y ~ 53. It is worth mentioning that 
shallow water flow has properties analogous to a two-dimensional gas flow 
with y = 2 (see Courant & Friedrichs 1948). 
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SUMMARY 

The decay of small perturbations on a plane shock wave 

propagating along a two-dimensional channel into a fluid at rest 
is investigated mathematically. ‘The perturbations arise from 
small departures of the walls from uniform parallel shape or, 
physically, by placing small obstacles on the otherwise plane 
parallel walls. An expression for the pressure on a shock wave 
entering a uniformly, but slowly, diverging channel already 
exists (given by Chester 1953) as a deduction from the Lighthill 
(1949) linearized small disturbance theory of flow behind nearly 
plane shock waves. Using this result, an expression for the 
pressure distribution produced by the obstacles upon the shock 
wave is built up as an integral of Fourier type. From this, the 
shock shape, €, is deduced and the decay of the perturbations 
obtained from an expansion (valid after the disturbances have 
been reflected many times between the walls) for € in descending 
power of the distance, ¢, travelled by the shock wave. It is shown 
that the stability properties of the shock wave are qualitatively 
similar to those discussed in a previous paper (Freeman 1955); 
the perturbations dying out in an oscillatory manner like ¢-**. 
As before, a Mach number of maximum stability (1-15) exists, 
the disturbances to the shock wave decaying most rapidly at this 
Mach number. A modified, but more complicated, expansion for 
the perturbations, for use when the shock wave Mach number ts 
large, is given in § 4. 

In particular, the results are derived for the case of symmetrical 
‘roof top’ obstacles. ‘These predictions are compared with data 
obtained from experiments with similar obstacles on the walls 
of a shock tube. 

1. INTRODUCTION 
The ability of plane shock waves, propagating into a stationary fluid, 
to retain their shape when subjected to small disturbances was considered 
theoretically in a previous paper (Freeman 1955). ‘This property has been 
called ‘stability’. It was realized however that the model used to demonstrate 
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this, namely, a perturbation to the shock wave produced by moving a plane 
piston impulsively from rest at constant speed into the stationary fluid, 
would be difficult to reproduce in the laboratory. In order to overcome 
this disadvantage another model is here studied. In this model, a plane 
shock wave propagating down a uniform channel into a fluid at rest is 
disturbed by small protuberances on the channel wall. ‘Chis model bears 
a close resemblance to actual experimental conditions when the shock is 
produced, as is usually the case, in a shock tube, and, in fact, such experiments 
have been performed by Mr K. C. Lapworth in the Manchester University 
shock tube. It will be shown that the shock wave qualitatively exhibits 
similar properties to those described in the preceding paper (Freeman 
1955). 

Chester (1953) has already considered the problem of a shock propagating 
along a two-dimensional channel which begins to diverge slowly, and has 
extended this result to obtain the change in strength of a shock propagating 
along a channel which changes from one uniform cross-section to another. 
The solution given below will be a further development of this work obtained, 
essentially, by considering higher order terms in the expansion given by 
Chester. In its turn, the Chester solution is a development of the problem, 
investigated by Lighthill (1949), of a shock wave propagating along a 
plane wall which suddenly changes in direction by a small angle 6. ‘Uhe 
Huid is assumed inviscid except in a thin region comprising the shock wave 
across which the variables of the flow change discontinuously in a manner 
determined by the Rankine-Hugoniot relations. ‘Vhe disturbances produced 
at the corner are assumed to be small enough for the equations of the flow 
behind the shock wave to be linearized. Vhe shock wave is disturbed from 
its plane form within a region cut off by a cylindrical wavefront, which 
originating at the corner, expands with the speed of sound behind the 
undisturbed shock wave and is translated bodily with the velocity ot the 
fluid behind the undisturbed shock wave (see figure 1). When the corner 
is concave to the flow, this becomes the familiar Mach reflection. Under 
the conditions stated above, the pressure satisties the wave equation with 
constant sound speed and the solution can be obtained in terms of the 
‘conefield’ variables by the use of the Busemann transformation and 
conformal mapping techniques. Chester has shown that the presence of 
another wall can be accounted for by considering a system of images in the 
wall. The solution is rather complicated, being a sum of Lighthill solutions. 
After the disturbances produced have undergone a sufficiently large number 
of reflections at the walls, however, an asymptotic expansion for the pressure 
on the shock wave can be obtained in descending powers of the time. ‘The 
first term of this expansion is used by Chester in his solution for a change 
in cross-section of the channel. Further higher order terms, which decay 
with time, are neglected. It is these terms that we shall consider here. 
In the case when the channel is uniform along its length, except over a 
finite range, after which it returns to its original cross-section, the decaying 


terms alone remain. 
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We, therefore, consider the propagation of a plane shock wave along 
a uniform two-dimensional channel with walls a fixed distance apart except 
over a finite length where the channel width varies slowly in some prescribed 
manner. Following Chester, we will concern ourselves with the behaviour 
at large time, that is, after many reflections have taken place at the channel 
walls. It will also be assumed that sufficient time has elapsed for the shock 
to have travelled a distance large compared with the length of the non-uniform 
section of the channel. 


Figure 1. The configuration in the Lighthill problem. 
S, shock; IV, wall. (a) Subsonic; (4) Supersonic. 


‘he perturbations to the flow arising from each element of the disturbing 
wall are propagated, as pointed out in the previous paper (Freeman 1955), 
at the speed of sound behind the undisturbed shock as cylindrical waves 
convected with the fluid. The velocity of the disturbances along the shock 
wave is the velocity component of these wavefronts along the shock. ‘l’o 
the linearized approximation, the velocity of the wavefront of the cylindrical 
wave produced by a particular element of wall is a, \/{1—([U—v]/a,)*} 
(= a, sind, say) along the shock, where a, is the velocity of sound behind 
the undisturbed shock and U—uw is the velocity of the undisturbed shock 
relative to the fluid behind it. As these waves are cylindrical in character, 
we expect a singularity like y {1 —(r/a,t)?} where r is the distance from 
the centre of the disturbance produced by the particular element under 
consideration, and ¢ is the time measured from the instant at which the 
shock strikes it. But, as has been shown previously (Freeman 1955), the 
strength of this wavefront is forced to vanish at the shock by the form of 
the boundary condition there, and the singularity is of higher order. The 
etfect of the non-uniformity of the channel wall is then the sum of all these 
cylindrical waves and is obtained as a Fourier transform of the individual 
solutions. At large time, the behaviour of this integral is dominated by the 
singularity at the cylindrical wavefront. ‘Thus, the shock is perturbed on 
passing the non-uniformity by disturbances which travel along it with a 
speed a,sinw. Due to the singularity at the wavefront, the disturbances 
deeay in an oscillatory manner like 3? where ¢ is the distance travelled 
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past the obstacle. ‘This rate of decay varies with shock Mach number, 
decreasing rapidly both for strong and very weak shocks. Between these 
extremes there lies a Mach number of maximum stability where the 
perturbations to the shock wave decay most rapidly with the distance 
travelled by the shock. ‘This has been calculated as MW = 1-15 approximately, 
which is very close to the value 1-14 given for the first model (Freeman 1955). 

An interesting phenomenon occurs when we consider the particular 
case of symmetrical ‘roof top’ obstacles placed on the upper and lower 
walls of the channel (see figure 2). ‘Then, together with the above variation 
in Mach number, there is a variation with the lengt! of the ‘roof top’ due 
to cancellation or reinforcement in varying degrees of the disturbances 
produced at the corners and vertices. A similar effect occurs in linearized 
aerofoil theory of supersonic flow for the so-called ‘Busemann Biplane’. 
When the shock wave strikes the ‘roof tops’, disturbances propagate 
outwards along it and are reflected at the opposite walls. At the ends of 
the ‘roof top’, the corners being concave to the flow, the disturbances will 


2b 


Figure 2. The ‘roof top obstacles in position in the channel. 


be compressive ones, whereas at the vertex, the corner being convex to the 
flow, the disturbance will be expansive. ‘The latter disturbance, according 
to linear theory, is opposite in sign and twice the magnitude of either of 
the former. On this account, should the length of the ‘roof top’ be such 
that the time taken for the shock to travel between corner and vertex is an 
integral number of times that for the disturbances on the shock to travel 
across the channel, cancellation of these disturbances will occur and the 
decay will be faster than indicated above. In a similar manner, maximum 
disturbance is produced when the time taken tor the shock to move from 
corner to corner is an odd integral multiple of that for the disturbances 
to cross the channel. ‘These conditions may be written Ud/a,/sinys = 2n 
and Ub a,/sind = 2n+1, where n is an integer, U is the velocity of the 
shock wave, 2/ is the length of the ‘ roof top’ and 24 the width of the channel. 

In $5, the results obtained by Lapworth (1956) from experiments on 
‘roof top’ obstacles situated on the walls of a shock tube are briefly discussed. 
These results indicate that the theory developed in this paper gives a good 
qualitative picture of actual conditions. The oscillatory nature of the 
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perturbations and their rate of decay are in good agreement although the 
actual magnitude of the perturbation is somewhat less than predicted. 
This discrepancy might be expected however since the piesent theory is 
only a linear approximation to the true state of affairs. 

In conclusion, a comparison of the present model with the one discussed 
previously (Freeman 1955) makes it clear that the behaviour of the shock 
wave is to some extent independent of conditions prevailing in the flow 
field behind. A theory making this assumption has been put forward by 
Whitham (1957) to take into account non-linear effects. 


2. EXTENSION OF CHESTER SOLUTION 

‘Lhe theory of Lighthill (1949) has shown that the flow behind a plane 
shock moving along a plane wall which changes in direction by a small 
gagle 6 is disturbed within a circle expanding with the velocity of sound 
behind the shock and whose centre moves with the velocity of the fluid 
behind the shock originally being situated at the corner (see figure 2). 
These velocities are, to this approximation, constant and equal to their 
value behind the undisturbed shock wave. In the case of supersonic flow 
behind the shock, there is an additional region of disturbance bounded by 
the tangent from the corner to this circle (figure 2(b)). The shock is 
disturbed from its plane form in the region cut off by this circle. 

The pressure within the circle is given by 


(Pi — Po) p(X, Y), (2.1) 

where 

Oop C[D(z,— 1+ y?)— 1]sec 

with cos’ = [(.W*+5)(7M*—1)]'*. The constants «, 8, y, C and D are 
functions of M detined by Lighthill (1949). .W is the Mach number of the 
shock; dp, Po, Py and a,, Py, p; are the sound speed, pressure and density 
in front and behind the shock respectively. ‘Vhe z,(= x, +7y,)-plane is 
related to the physical (x, v)-plane with origin at the corner and x-axis in 
the direction of propagation by 


(2.2) 


1 pen" +y) ] ) 2 if — ell 
{ pe’ — f perv) f | (2.3) 
and v—-g,t .v 2pe'” 
= Yigy = p 
a,t a,t 1+p? 


where q, is the velocity of the fluid behind the undisturbed shock wave. 
The disturbed region in the (x, y)-plane is mapped into the upper half plane. 
The disturbed shock, which is given by x = Ut ar VY = cosy, becomes 
1+(¥ sin)? 

1—(Y/sin ys)?’ 


with Y = sing. 


= ¥, = 0 (2.4) 
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The point x, = 1, y, =0 is the intersection of the shock and wall and 


v, = “&, vy, = the intersection with the circular wavefront. From (2.2), 
cp C[D(x, — 1+ y?)— 1](a+ B)sece (2.5 
dx, (x, + 72x, 1+ -14 22) 


on the shock. 

Chester (1953) has shown that the problem of a shock moving along a 
channel of width 24 which begins to diverge uniformly at a small angle 26 
can be solved in terms of this solution. ‘The effect of enclosing the shock 
between two walls is to reflect the waves produced at the corners repeatedly 
between the two walls. Since this wave motion is already of order 4, 
however, the effect of the divergence of the walls will be a second order 
effect in 6 in the reflections, and hence these reflections may be considered 
to take place at plane walls. Suitable images can be placed outside the walls 
to produce the required reflections. ‘The number of these images that 
influence the shock will grow with time. The pressure in the case of a 
channel can therefore be written 

(Pi — Po) P(x, 
with 


v—qg,t (2n+1)b-y v-g,t (2n+1)b+y 
a, t a, t at f 


where x is measured along the channel and y from the centre of the channel, 
the summation being taken over all the images that influence the shock. 
By using Dirichlet’s summation formula, Chester writes this, in a torm 
more convenient for asymptotic estimation, as 


P = a, t)+ Lp(X, a, t)- 


+ | (—1)"p(X, (2.7) 
~ 4) ayt 
where p(X, }) is given by (2.2). 
If we now consider pressure variations on the shock alone, then 


P=] ( , ory) 
= } “OS Ws, us, = + 


a,t 
| (—1)"p(cosw, jentiian t¥ gy 


1)" p(cos us, Y jer tien t¥ y) b (2.3) 
b+yayl 

It will be suthcient to consider the behaviour of this function for large time. 

The behaviour of the first two terms will simply be determined by 

expanding p as a function of }. ‘The integrals within the sum will however 

be dominated by the singularities of p in the range of integration. It will 

be shown as in the previous paper (Freeman 1955) that the singularity at 


the shock-wavefront intersection 1s dominant. 
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Let us now consider the function p(cos¥%, Y) in more detail. Near 
the wall, when Y = 0, x, = 0, (2.5) together with (2.4) states that 
p = const. + O(Y?). 
Thus cp/cY = 0 at Y = 0 and also 
p(cos&, (b+ y)/a, t) = const. + O((b + y)?/a? (2.9) 
At the shock—wavetront intersection Y = sing, or x, = 2%, then, from (2.5), 
op CD(x + B)sec 


Ox, 


(2.10) 


Consider now a general term of the sum in (2.8). For m 4 0 we have 


J, = >> (—1)"p(cos&, dy 
2b (b¥ wlayt 
—1)"¢ ‘sin Cp. 
| J ay t 


using (2.9) and (2.5), where 4 is a constant. Also, since ep/eY = O(Y) 
as Y +0, it follows that 
( rsiby op 


n 


t¥/b (5). (2.12) 


where B,, is a constant and the integration is now along the shock from the 
wall to the wavefront. It follows also that the above integral is dominated 
by the singularity at the upper limit. The integrand is 


shock 


~ CD(x+ 8)sec cosec — ¥ cosec (2.13) 
near Y = sind, and so the integral may be written 
CD(xz+ 8)sec us cosec | (1 — cosec ys)! + 


Using the formula for the asymptotic expansion of a Fourier integral (see, 
for example, Freeman 1955), we obtain 


nza,t | 
ina 
exp| (a, tsin’ + | (2.14) 
where (n>), 
= +] (n < U), 
for ¢ large. Substituting in (2.8), we have 
C D(z + 8)sec 32 —j)" 
D(x d > ( 
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where >’ denotes that there is no term m = 0 in the summation, and B is 
a constant. ‘The integral has been calculated by Chester, who denotes it 
by A(U/a,), and hence 


KUt RB + B)sec ys bh 32 
a,tsing) 
2 cos( = ) (2.16) 


Now, let us consider a channel of more general cross-section defined by 
vy =hb+f(x), v = —b—f(x), where f(x) is non-zero for some range of x and 
f(x) is small. In fact, we will consider the case of f(x) non-zero only in 
a finite range, that is, the effects of change in shock strength due to changes 
in width of channel, which have already been studied by Chester, will not 
be considered. ‘The solution of the problem in which the channel begins 
to diverge uniformly can now be generalized to this new channel shape by 


| considering a change of slope df'(y) at time t = 7/U and distance x = 7. 
q Thus the value of P for the new problem will be 


the integral being taken over the range in which f is non-zero. For 
sufficiently large distances away from the obstacle in the channel, the 
powers of t—(y LU’) become ¢ and thus, 

CD(x- 


5 
a, tsinis 


| \ 


¢(-1)" nza, sings 
| n° * cos h (: U 4 


Ty l 
x cos df'() + (3): (2.18) 


Now if € is the perturbation of the shock wave from its plane form, then the 
perturbation to the pressure is 2.Mayp)(c& ct), to a linear approximation, 
from the shock equation (y = 1-4). Hence the perturbation to the shock 
shape may be written 


3CD(a+ B) (2 Py ( h 
h 573M cos sin ys\ a, Ay py] t sin ds 


2 ¢(-1)" Pnza,sind’ 7 ,., 
> S 1 — + (2.15 


Vhis result indicates, as has been shown previously (Freeman 1955), that 
the perturbations to the shock wave decay at a rate proportional to the 
inverse three-halves power of the time, or, alternatively, of the distance 
i travelled by the shock wave past the obstacle. The variation of the coefficient 
ze in (2.19) with Mach number will be discussed in more detail in § 3. 
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In passing, it might be observed that the modification to equation (2.19) 
tor the case of a channel which increases in width by a small amount 2c is 
obtained by adding a term 

— po) K, Ut 
5a, US 
to the right-hand side. 


3. "THE ‘ ROOF TOP ’ DISTURBANCE 
In order to obtain more detailed knowledge of this phenomenon, we 
will now study a particular form of the function f(*)—the so-called * root 


top’ (see figure 2). Let 


(-1<9 <0), | 
| 
| 
| 


(3.1) 
Then, using (2.17), we obtain 
P, 2P (ty) P(1 (3.2) 


The total disturbance is made up of a sum, therefore, of the individual 
disturbances at the corners and vertex of the ‘roof top’, due regard being 
taken of their origin in time. Hence, using (2.19) and (3.2), 


bh Na tsing) ~ 
| 28(5 (a, tsinds +») + (a, tsinyds 
7 
U) +y)) (asingir+( U)} -»)- 
(asin uit—(1 (asin -») |. (3.3) 
where Eiz)= > (sin ns — cos ns) 
n=1 
and Q, = 3CD(z+ BYP, — Po) 


107%. cos sin ba, dy po 


If ¢ = Ut, the distance travelled by the shock, then 


— 7a, Ising _ ple’. za, 
where F(z) = 2E(z) (3.4) 


U 32 
and Q 1 


es 

at 


406 N. C. Freeman 


Figure 3. The function F(z). 
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Figure 4. The coefficient {2 as a function of M. 
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‘The function E(s) was computed by making use of the relation 
| 0 1 


| 


The function F(z) is plotted in figure 3 for various values of /a, sinw/ Ub. 
F(z) is periodic with period 27. For /a,sinys/Ub = 2n, where n is an 
integer, the function is identically zero. ‘This condition may be written 
LU = (26/a, sinys)n, which states that the time taken for the shock to travel 
between a corner and the vertex of the ‘roof top’ is an integral multiple 
of the time taken for a disturbance to travel along the shock across the 
channel. If this is so, the initial term of our expansion will vanish and 
higher order terms will dominate, the shock perturbations decaying more 
rapidly. ‘The function F is a maximum for values of (a, /sin & Ub) at which 
(//U) = (ba, sin&)(2n+1). These values correspond to configurations 
in which the disturbances produced at the vertex and corners reinforce 
each other. ‘The coefficient Q is shown in figure 4. The similarity between 
this figure and figure 5 of the earlier paper (Freeman 1955) is immediately 
obvious. The minimum at M = 1-15 is very close to the one shown there. 
The perturbations to shock waves of this Mach number decay most rapidly 
with the distance travelled by the shock wave. 


4. "THE VARIATION AT HIGH MACH NUMBER 
‘The behaviour of the various constants for large Mach number M may 
summarized as follows: 
sind~y(%), B~ \ (3), 
The correct approximation to equation (2.5) for large M and large x, 
is, therefore, 


b 


fa) 


Ox, a J + [(x, 

since x% is no longer small. It will be observed that as x « then 
(cp/cx,) ~~ x °? and the form of the singularity that dominates the integral 
(2.13) has changed. When « is large but finite, we must consider a form of 


expansion in which the pole x, = —? lies close to the branch point. ‘The 
method used is due to Clemmow (1950). From (4.2), 

(4.3) 


near s = 0, where = and O = CD(xz+8)secy Whence, 
using (2.4), 

Cp 

sin a-*)’ 


(4.4) 


and thus for large WV, 


Cp () 
éY ~ sing{1—(¥ sind)? 
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as opposed to op 

sinds)}!? 
for smaller /. ‘This in turn would give a slower decay, like t-!?, to the 
shock perturbations for large M. For large M, therefore, we require an 
expression for 


(4.5) 
Jo +a? 

in place of the expansion of (2.13) since, as before, the main contribution 

will come from near s = 0. ‘Thus, as in §2, the pressure due to a slightly 

diverging channel is 


( ) | einem, tsin dz inaa,tsiny/b O (4.6) 


The integral enclosed in the square brackets may be written 


na, t sing xb 


and b \'2 / n za,tsingys 
fi fora < 0, 
na,tsinds 


where * denotes the complex conjugate and 


J 
As in (2.19), the shock shape is then given, for a channel varying in width, 
by 
\ 


(¢ bh)! h 
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where \ is given in table 1. And, in a similar manner to §3, the solution 
for a ‘roof top’ obstacle can be obtained. 


| M \ | 
| 
1-5 0-102 
20 | 0-128 | 
3-0 0-168 | 
x 0-234 


Table 1. 
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5. EXPERIMENTAL RESULTS 

The model discussed in the previous section bears a close resemblance 
to actual experimental conditions in a rectangular shock tube when ‘ roof 
top’ protuberances are placed on two opposite walls. Such experiments 
have been undertaken by Mr K. C. Lapworth (1956) in the shock tube at 
the Mechanics of Fluids Department, Manchester University. The original 
intention was to make a comprehensive study of the results obtained from 
the preceding theory. However, the results obtained so far fall somewhat 
short of this ideal due to the inability to reproduce the theoretical model 
as closely as would be desired. 

In theory, the perturbations to the shock wave are assumed infinitesimally 
small—a consequence of the linearization of the problem. In practice, it is 
necessary to produce disturbances which are large enough to be measured. 
Since the perturbations are decaying all the time, this will always limit the 
distance downstream at which one can observe the perturbations to the 
shock wave. As the theoretical result is given by an asymptotic expansion 
in inverse powers of the distance, this can seriously limit any comparison 
between theory and experiment by not allowing time for the regime in 
which this expansion is valid to be set up. In fact, however, the main 
features of the asymptotic form of decay would seem to be apparent earlier 
than one might expect theoretically, at least for Mach numbers not too small. 

Also, in the theoretical considerations viscous effects are neglected entirely 
except in so far as they occur in the idealized form of the shock wave. 
In practice, viscous effects are important especially around the corners of 
the ‘roof tops’, producing separation and the formation of a wake from 
the vertex. As any viscous effects produced at the wedges will take time to 
develop after the shock has passed, it will take time for such disturbances 
to grow large enough to influence the shock. It is thought that it is this 
which puts an even greater restriction on the distance downstream at which 
measurements of the perturbations under theoretical flow conditions can 
be made. For, after a time, the perturbations start to increase due probably 
to the wake formation influencing the shock wave. 

‘The experiments were made at three shock Mach numbers, 1-60, 1-41 
and 1-16. ‘The shock wave was photographed using Schlieren techniques 
at various distances after passing the ‘roof top’ system. As it was difficult 
to decide where the undisturbed position of the shock wave would be and 
hence to measure € directly, a measure of the perturbation which was 
independent of this was used. ‘This will be called ‘the total perturbation’ 
and denoted by €,. &, is the sum, without regard to sign, of the maximum 
and minimum values of € along the shock wave at any particular time. ‘The 
Mach numbers 1-41 and 1-60 give disturbances to the shock wave which 
are near the maximum for the ‘roof top’ system cgnsidered, and hence, 
for these, it is sufficiently accurate to approximate the function F(z) by a 
sinusoidal function of the appropriate amplitude and period. ‘The total 
perturbation €, can then be written in the form 


5 
= sin(mx +p), (5.1) 
F.M, 
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from equation (3.4), where the values of n, G, m and p are functions of 
Mach number and are given in table 2. 

For the Mach number 1-16, it was found difficult to compare the values 
with theory in any way. ‘This may be due to the fact that the asymptotic 
behaviour will take a longer time to establish itself for lower Mach numbers 
and those results will therefore not be discussed further. For the two 
Mach numbers 1-41 and 1-60, Lapworth (1956) has attempted to fit the 
experimental results for a formula of the type (5.1). Estimates of m and p 
are obtained from the zeros of the function €, and then G and » are obtained 
by applying a least squares fit. In order to decide at what point the 
perturbations began to be influenced by the disturbances in the flow 
produced by separation at the roof tops, the data were analysed using 
varying numbers of results dispensing one at a time with the later results 
until a fairly consistent set of values on n and G were obtained. An average 
of these values is given in table 3. 


| | 
M | n | m p | G | 
| 
| 
| | | 3-39 3-91 | 0-157 | 
| 1:60 455 3-44 3:91 | 0-273 | 
| | 
Table 2. 
- | | | | 
M | | m G 
| 145 3-31 | 0-066 
1-60 | 165 | 3-31 19 | 0-090 | 
Pable 3 


A comparison of tables 2 and 3 shows that even with the difficulties 
encountered, the decay of the perturbations is in fair qualitative agreement 
with theoretical predictions. ‘he actual perturbations are two to three 
times sinaller than those predicted but this might be expected from a linear 
theory. Non-linear effects would interfere with the reinforcement or 
cancellation of the sets of disturbances from the two roof tops and also 
cause the disturbances to interact among themselves. The compressive 
disturbances would tend to concentrate themselves into shock waves 
whereas the expansive ones would spread out and interact with these waves. 
A non-linear theory, such as that developed by Whitham (1957), would 
be required to describe these effects. The resulting perturbations would 
thus be reduced by a further cancelling of the individual disturbances. 
The poor agreement of the phase difference p may again be due to the 
non-asymptotic nature of the flow. 
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In view of the reason for the breakdown of the theoretical flow conditions 
in practice, it would seem difficult to modify this model in any simple way 
to overcome this defect. ‘To make any more detailed study of the phenomena, 
it would seem that a different shaped obstacle not so prone to separation 
difficulties would have to be used. 


I would like to express my sincere thanks to Professor M. J. Lighthill 
for his help and encouragement in the preparation of this paper, and to 
members of the Departments of Mathematics and Fluid Mechanics with 
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Lapworth, Dr D. Blackburn and Dr F. Goldsworthy. 

I am grateful to the Department of Scientific and Industrial Research 
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REVIEWS 


Relaxation Methods in Theoretical Physics, Volume II, by R. \V. 
SOUTHWELL. Oxford: Clarendon Press, 1956. 274 pp. 55s. 


It is net often that a scientific treatise is at once a comprehensive 
account of the subject and a record of much of the professional career of 
the author. But Sir Richard Southwell has made the subject of relaxational 
metheds peculiarly his own, and the work reviewed here fits the dual 
description. 

Anticipations or hints by early writers usually make it impossible to set 
precisely the date of invention of anidea. In 1924 Hardy Cross introduced 
one of the major concepts of the relaxational approach in conjunction with 
his ‘moment distribution’ method for the calculation of the forces in 
statically indeterminate structures. ‘Thom’s device for solving Laplace's 
equation (1928) ditfers from what we now call relaxation only by a small 
yet important technicality. Be this as it may, the concept and use of the 
relaxational technique as a systematic device for obtaining numerical 
solutions to physical problems of all descriptions derives from the investiga- 
tions of Southwell beginning in 1935. Since that time he has been the 
high priest of relaxation: a faithful worker, an enthusiastic proponent 
of its virtues, and, in a small way, the founder of a school. ‘The major 
part of the literature on the subject is by Southwell and by younger men 
who entered the field under his tutelage. 

This, the second volume of Relaxation Methods in Theoretical Phystcs, 
appears ten years after publication of the first. ‘he two together make 
a whole, and they have been so treated for the purpose of this review. 
In the preface to the first volume, Southwell conceived the work to be 
“of less concern to engineering than to theoretical physics”’, in 
contradistinction to his 1940 volume on Relaxation Methods in Engineering 
Science. Either Southwell has misjudged the situation, or in the meantime 
the scope of engineering has grown mightily, for most of the problems he 
deals with are of first-rate significance in engineering. The stated aims 
of the two volumes are to set out the underlying theory in a coherent 
and detailed manner, to describe the |. ious techniques of computation 
proved by experience, and to convey the rather different atitude conferred 
by the approach to a physical problem through relaxation. All this 
Southwell has done in a lucid if rather formal manner. 

Both volumes are concerned with physical problems of the classical 
continuum type. With one minor digression, the range of problems 
is confined to those having two independent variables: either two space 
coordinates, or one space and one time coordinate. The first volume 
begins with a statement of certain aspects of the calculus of finite differences, 
as a preliminary to the reformulation of partial differential equations as 
algebraic difference equations. About half of this volume is concerned 
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with the numerical solution of Laplace’s equation in the plane, including 
conformal mapping by numerical methods, and with the quasi-plane- 
potential problem. ‘The first volume closes with * free surface ’ problems 
involving boundaries or interfaces not initially known. The second 
volume continues with more difficult problems, many of which were still 
classified when the first volume appeared. It deals principally with the 
biharmonic equation; with eigenvalue problems; with non-linear problems 
of large elastic distortions, viscous fluid motions, and plastic straining; 
and it concludes with a brief resumé of recent innovations in relaxation. 

When a difterential equation purporting to describe a_ physical 
phenomenon does not yield to analysis, relaxation, with its entirely different 
outlook, may be used to advantage. The relaxer does not think in terms 
of independent and dependent variables as such. He fastens attention 
on the numerical values of the wanted quantity at the nodal points of a 
finite net in the plane (or the space) of the independent variables. 
Relaxation comprises a systematic method for determining these numerical 
values so that they satisfy the finite-difference version of the differential 
equation. Once the problem has been formulated in terms of relaxation 
operations indicating the steps required to liquidate the residues (or errors) 
of the finite-difference equation, its numerical solution requires only 
simple arithmetical calculations. 

Southwell quite rightly brings in much physical argument, for this 
helps the relaxer to get on with his job. But the reviewer is not the first 
to be irritated and distracted by Southwell’s insistence on thinking of so 
many problems in terms of the analogy of the force-displacement system 
in a net made up of strings. ‘This does in fact permit one to revive in 
a rather strained form the original meaning of the term ‘ relaxation’, 
that is, the systematic relaxation of artificial constraints at the joints of 
a redundant structure, until, all such artificial constraints having been 
relaxed, the solution sought for is at hand. But while the string-net 
analogue may somehow illuminate for Southwell problems in, say, potential 
fluid motion, it interferes with clarity of presentation for others. 

Relaxation has a brute force quality which makes it possible, at least 
in principle, to solve numerically any problem that can be quantitatively 
formulated in terms of algebraic, differential or integral equations. The 
wielder of the pencil need know no mathematics beyond arithmetic. 
Nevertheless mathematical skill and physical insight are assets which 
usually greatly reduce the time required to effect a solution. Above all, 
the calculator should have the qualities of an accountant—patience, 
neatness, and an empathy for numbers. While the arithmetical operations 
are simple enough, the whole process proceeds rapidly only after a certain 
amount of experience has been gained. On new and exceptionally difficult 
problems, the circumstances are often such that various bits of analysis 
must be invented ad hoc as the calculations proceed, and the formulator 
of the problem may easily find himself in the position of doing the tedious 
arithmetic. After the few weeks or months required to complete his 
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problem, he has become an expert relaxer, but unless he is one of those for 
whom relaxation possesss a special fascination, he may thenceforth avoid 
it like the plague. 

The present position of relaxation methods so far as the differential 
equations of mathematical physics are concerned is that it has indeed 
demonstrated its ability to solve problems before which more orthodox 
mathematical methods fail. From a practical point of view, these may be 
classified in two categories. In one class, e.g., Laplace’s equation for 
plane, irrotational motion of an incompressible fluid, the calculations are 
so straightforward and progress so rapidly toward a solution that the 
method is perfectly feasible for the engineering office or for almost any 
circumstance where, say, a day’s labour is not excessive. ‘The other class, 
e.g., the motion of compressible gases at high speed, or the viscous motion 
of fluids when inertia cannot be ignored, contains problems requiring the 
labour of weeks or even months, with the frequent application of analytical 
thinking. Such problems demand a determined and motivated individual 
who can spend the necessary time without distraction. One need hardly 
say that such is not likely to be the case in an engineering office. But 
the solutions of even a few such problems may settle certain crucial questions 
arising in a research effort, and each problem solved may exemplify some 
phenomenon in such manner as to enlarge greatly our understanding if 
not our capacity to deal quantitatively with it. 

An enumeration of titles from the ‘Index of Problems Solved” 
delineates the range of problems in fluid mechanics that have yielded to 
relaxational methods : 

Conformal transformation of an area into a circle; 

Conformal representation of the region external to an airscrew 
section into the region external to a circle; 

Oil pressure-distribution in a plane slider bearing; 

Oil pressure- and temperature-distributions in a_partial-journal 
bearing; 

Flow of gas through a convergent-divergent nozzle; 

Percolation through an earth wall; 

Laminar flow through a two-dimensional Borda mouthpiece ; 

Orifice plate in a circular tube; 

Free jet falling under gravity (* waterfall ’); 

Slow motion of viscous fluid through a two-dimensional nozzle; 

Only a fraction of the book is concerned with problems of fluid 
mechanics, however. Major sections deal with elasticity, plasticity, and 
vibrations, and a few problems in heat conduction and electricity are 
discussed. 

The final chapter of Volume II, apparently written after the main 
body went to press, gives brief accounts of two recent ideas, both due to 
D. N. de G. Allen. ‘The first employs the device of isometric representation 
on a single sheet of paper to effect a three-dimensional relaxation solution, 
a feat which would strain the sanity of any one attempting it on multiple 
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leaves of paper. The second converts a ‘ marching’ type of problem 
to one of ‘ jury ’ type, thereby making the problem accessible to relaxation 
methods. ‘The trick here is to effect a substitution of variables which yields 
an equation of higher order; since the actual number of physical boundary 
equations is unaltered, one or more boundary conditions in the substituted 
variable becomes free both in its value and in its location, and the problem 
expressed in the substituted variable may thus be framed in the jury form. 
By this device the parabolic differential equations of unsteady heat conduc- 
tion, for instance, may be solved by relaxation. 

Although Southwell’s single-minded enthusiasm might lead one 
to believe otherwise, there are few circumstances in which a relaxation 
solution would be preferred to an analytical solution. Questions of 
computation time aside, the analytical solution possesses two decisive 
advantages. A single analytical solution in closed form actually comprises 
many solutions by reason of the parametric constants which appear. 
And, of great importance to the applied physicist and engineer, it exhibits 
in condensed form certain patterns of physical behaviour which no mere 
collection of sheets of relaxation solutions, each with its own array of 
numbers, can possibly evoke. 

Southwell’s presentation of the subject incorporates a balance between 
the physical, the mathematical, and the computational sides. Much 
of the subject matter is discussed in reference to specific physical problems, 
and many complete numerical solutions are displayed elegantly and 
clearly. However, one looks in vain for any critical assessment of the 
physical models employed, or for meaningful interpretations of the results. 
Southwell considers the problem ‘ solved’ when the numerical solution 
is completed. ‘This should be taken as an observation rather than a 
criticism, for the book is really concerned with relaxation methods and not 
with physical problems. The mathematical aspects are given their full 
share of attention, with the thought that this would be one of the more 
useful items for those who might wish to apply relaxation to an as yet 
unsolved problem. ‘lhe various labour-saving tricks of computation and 
schemes for economically organizing the work are discussed to the point 
where further elaboration would not be as useful as actual practice on the 
part of the reader. 

It is disappointing that, in a book concerned with numerical methods 
of solution, there is no discussion of the relationship of relaxation methods 
to the recent flourishing progress in high-speed machine computation. 
Relaxation calculations for differential equations with two or three 
independent variables pose difficulties for a high-speed machine because 
information must be remembered at a great many net points. Several 
questions come to mind. Having regard for the, special characteristics 
of digital computers, how can the relaxation procedure be adapted to use 
best the machine’s limited memory but great speed? Is it feasible or 
desirable to use many fewer net points in conjunction with higher-order 
terms in the Taylor-series expansion about a central point? Is it better 
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for the machine to search for the largest residual or to relax the residuals 
at the net points in some fixed geometrical sequence? The reviewer 
would have welcomed reflections along these lines. 

Finally, brief mention may be made of two other books on relaxation 
methods*. ‘These are both much shorter books dealing with the whole 
of the subject in the sense that algebraic equations and ordinary differential 
equations, as well as partial differential equations and eigenvalue problems, 
are treated. Shaw’s book omits all discussion of physical problems, and 
is focussed on the mathematical theory underlying relaxation and on the 
computational procedures. Allen’s book is more akin to Southwell’s 
in its description and illustration of the subject through the medium of 
physical problems. 

For the novice to the subject, Allen’s book is, in the opinion of this 
reviewer, the easiest and best organized introduction to the subject. 
Southwell’s book is the definitive treatise on the subject, giving extensive 
details of the physical, mathematical, and computational sides, and 
constituting a compendium of the significant types of problems accessible 
to relaxation. And Shaw’s book seems best suited as an introduction for 
those interested only in the mathematical and computational aspects of 


relaxation. 
A. H. SHAPIRO 


Proceedings of the 8th International Congress on Theoretical and 
Applied Mechanics (Istanbul, 1952). Published by the Faculty 
of Science of the University of Istanbul. Vol. 1, 1953; 529 pp. 
Vol. 2, 1955; 203 pp. £3 15s. Od. for the two volumes. 


These two volumes providing a record of the 8th International 
Congress for Applied Mechanics, held at Istanbul in August 1952, are 
now available in Britain and U.S.A. The first volume contains the list 
of members of the congress, the scientific programme, the addresses read 
at the opening and closing sessions, and summaries (a page or two in 
length) of the 336 papers read at the congress. ‘The second volume 
contains extended accounts of most of the seven general lectures and ten 
sectional addresses. 


* An Introduction to Relaxation Methods, by F.S. Shaw. New York: Dover, 1953. 
Relaxation Methods, by D. N. de G. Allen. New York: McGraw-Hill, 1954, 
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